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PREFACE 


This document contains a collection of 21 papers presented at the AIAA/ASME/ 
ASCE/AHS/ASC/ 31st Structures, Structural Dynamics and Materials Conference 
held in Long Beach, California, April 2-4, 1990. The conference had a total of 
284 papers, including 263 full-length papers and 21 short presentations in the 
work-in-progress sessions. All of the papers appearing in this document were 
presented in the two work-in-progress sessions. Most of the full-length papers 
are contained in the conference proceedings published by AIAA. 

The fields covered by the conference are rapidly changing, and if new results 
and anticipated future directions are to have maximum impact and use, it is 
imperative that they reach workers in the field as soon as possible. This 
consideration led to the decision to publish these proceedings prior to the 
conference. Special thanks go to the members of the Research Information and 
Applications Division at NASA Langley Research Center for their cooperation 
in publishing this volume. 

The use of trademarks or manufacturers’ names does not constitute endorsement, 
either expressed or implied, by the National Aeronautics and Space Administration. 


Jean-Fran 5 ois M. Barthelemy 
Ahmed K. Noor 

Compilers 
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Test/Semi-Empirical Analysis 


Since 1975, extensive testing of carbon/ epoxy tape plates and 
stiffened panels has been performed (Reference l through 6) . 
Attempts were made to predict the crippling failure of stiffened 
panels, fabricated from C/Ep tape, using the non-linear option in 
the STAGS computer code (Reference 7) . However, no meaningful 
results were acquired. Therefore, a semi-empirical crippling 
method was developed. 

To date, a semi-empirical analysis method has not been developed 
for plates and stiffened panels manufactured from C/Ep fabric. 

The purpose of this work-in-progress is to present a semi- 
empirical analysis method developed to predict the buckling and 
crippling loads of carbon/epoxy fabric blade stiffened panels in 
compression. This is a hand analysis method comprised of well 
known, accepted techniques, logical engineering judgements, and 
experimental data that results in conservative solutions. In 
order to verify this method, a stiffened panel was fabricated and 
tested. Both the test and analysis results are presented. 
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Bucklinc 


Test. Specimen 


This figure shows the test panel configuration. It consists of a 
skin with three blade stiffeners. The blade stiffeners contain 
flanges which were cocured to the skin. The entire panel was 
made from Hercules AS4/3501-5A carbon/ epoxy fabric except for the 
C/Ep tow used at the flange/blade intersection. This C/Ep tow 
provides structural integrity at the joint, including significant 
torsional stiffness provided at the blades. 

The blade stiffened panel was completely A and C-scanned and no 
defects were found. Prior to test, the panel was machined and 
assembled with potted aluminum end channels. The end surfaces 
were then ground parallel within .001 inch. 

The unloaded edges of the outer skin elements were supported by 
split rigid steel tubes to simulate simple support boundary 
conditions. This isolates the three stiffeners as though they 
were in a much wider stiffened panel. Thus, it was sufficient to 
analyze just the middle stiffener and apply this result to all 
three. The load carried by the skin adjacent to each split tube 
was justifiably neglected because it is such a small percentage 
of the total panel load. 



Aluminum Channel (Typ) 


Split Steel 
(Held on \ 
Not S 
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Detail A - Cross Section 


This is the cross section of the stiffener/ skin intersection. As 
mentioned above, the region adjacent to the blade, between the 
flange & skin, was filled with longitudinal carbon/epoxy tow. 

This juncture provides substantial support to the skin and the 
blades. However, the load carrying capability of the tow is 
neglected in the analysis. 

The panel elements were configured so that the skin buckled first 
and the blades buckled second. Thus, the flanges, which buckle 
last, support both the skin and the blades. 
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Typical Blade/Skin Intersection 


This is a photomicrograph of the manufactured stiffener/skin 
intersection. Good consolidation was achieved and structural 
integrity of this joint was expected. The curvature of the blade 
middle plies was inadvertent, but no reduction of boundary 
constraint was predicted. 



No-Edqe-Free Postbuckl inq Test 


In order to develop a semi-empirical stability analysis for 
carbon/epoxy fabric stiffened panels, empirical buckling and 
crippling curves for plates were generated. The plates tested 
were symmetric and balanced C/Ep fabric laminates. Each test 
plate was rectangular with clamped boundary conditions on the 
loaded edges (i.e., the short sides). Various b/t ratios were 
examined. 

Two unloaded edge boundary conditions were tested. The first, 
designated "no-edge-free" , was simply supported on both unloaded 
edges. The second, designated "one-edge-free" , was simply 
supported on one edge and free on the other. 

This is a typical no-edge-free plate test in compression. The 
unloaded edges are supported by steel v-blocks, simulating 
simple-support boundary conditions. The test specimen is in a 
postbuckled state. A full longitudinal wave can be seen. 
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PI 


Postbuckling failure of the no-edge-free compression test 
specimen is shown. This failure is referred to as "cripp 
The type of failure shown is typical for carbon/epoxy fab 
plates . 
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One-Edqe-Free Postbuckl inq Test 


This is a typical one-edge-free plate test in compression. One 
unloaded edge is supported by a steel v-block, simulating a 
simple-support boundary condition while the other unloaded edge 
is free. The test specimen is in a postbuckled state. One 
longitudinal half-wave can be seen. 
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One-Edqe-Free Crippling Test 


Postbuckl ing failure (or crippling) of the one-edge-free 
compression test specimen is shown. This type of failure is 
typical for carbon/epoxy fabric plates. 
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Typical lpad-Displacement Curve from Crippling Test 


This is a typical load-displacement curve for a plate compression 
test, either no-edge-free or one-edge-free. Displacement refers 

to the end-shortening of the test specimen. Buckling (P cr ) 
occurs at the bifurcation point of the linear curve. Crippling 

(P ) is the maximum postbuckling load that is reached prior to 
failure. 


Load 

(Kips) 
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No-Edae-Free Buckling Graph 


The no-edge-free buckling test data shown defines an empirical 
buckling curve for composites with similar layups. The ordinate 
is the ratio of the test buckling stress divided by the 

calculated classical buckling stress (F cr /Fcl^ * The abscissa is 
the width-to-thickness (b/t) ratio. The value for the classical 
buckling strength (F°£) can be obtained by using one of the 
following eguations. 

* Simply Supported Unloaded Edges 


F cr,u,0E 
cl, i,ss 


2jr, 

tb‘ 


£ ^ D 11 D 22^ 


+ D 


12 


2D 66> 


* Fixed Unloaded Edges 


cr , u , <pE 
cl , i, fx 



[4.6(D 11 D 22 ) 


+ 2 . 67 ( D^2 ) + 5 . 3 3 ( Dgg ) ] 


The classical buckling stress can be quite unconservative at low 
b/t ratios. However, the classical theory is accurate at b/t 
ratios greater than 50 . 


Fcr/Fcr 
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K hS i,?? e-ed ^ e free buckling test data shown defines an empirical 
SfiJS 5 urve f° r ? llnilar composite layups. The value for the 
us?nS d ?bI^ free ciassrcai buckling strength can be obtained by 
using the following equation. 2 


F cr , u, IE 
cl , i , ss 


t(L') 


where L' = 


eauJi^J ? n « J lxit y ^°® ffi cient of columns and is approximately 
equal to 3.6 for potted end columns in a test machine. 

This graph and its use is similar to that for no-edge-free 
composite plates. The discrepancy between classical and 
experimental buckling at low b/t ratios is the result of low 
transverse shear stiffness (Reference 8) . This effect is 
insignificant at large b/t ratios. errecc is 


Fcr/pcr 

cl 



b/t 
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Laminate Ultimate Compressive Strength 


• cu 

This figure shows the ultimate compressive strength (F ) for 
AS4/3501-5A fabric 0°, 45° composite laminates. This data was 

generated because F is required for the nondimensional 
empirical crippling curves which follow. 



AS4/3501-5A Fabric 
0°, 45° Laminates 
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The no-edge-crippling test data shown was used to define the 
approximate mean crippling graph. The ordinate is the crippling 

stress (F cc ) divided by the classical buckling stress (F^) , 

while the abscissa is the ultimate compressive stress (F cu ) 
divided by F cl . Thus, for plates with similar layups, where the 

value for F cu is known and F^J can be calculated, the predicted 
crippling stress may be obtained. 


Fcc/Fcr 
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Crippling Strength Predictions 


Armed with the empirical buckling and crippling curves, a step by 
step process can be used to calculate the crippling strength of 
the middle stiffener's blade and flanges. The classical buckling 

• C1T S t 

strains (e* ' ) for the blade and flange elements are .00309 

in/in and .00582 in/in, respectively. In this case, the blade 
buckles first and causes the flange element to buckle 
prematurely. Therefore, the minimum classical buckling strain is 
equal to .00309 in/in. 

The classical buckling strain is theoretical, not actual, and is 
referred to as a pseudo-strain. Using this strain (.00309 in/in) 

and the elastic modulus (E^'V) , the pseudo-buckling stresses 

(F* r, V) are calculated. The compression strengths ( F ^ U/ ^) are 

found from lamination theory or test results. The resulting 

pseudo-crippling stresses (F. '.) are then obtained from the one- 

/ 1 

edge-free empirical crippling curve. 

However, the empirical crippling curve was developed from testing 
plates with simply supported boundary conditions. The actual 
blade has a boundary condition better than simply supported but 
certainly not fixed. Consequently, the boundary condition was 
assumed to be equal to one-half the increase in fixity from 
simply supported to fully fixed. This correction factor (Ca) is 
only applied to the blade because the flange supports the blade 
until it buckles, but at that point, the blade cannot provide 
greater than simple-support to the flange. The crippling load of 
the middle stiffener is obtained by summing the stiffener element 

pseudo-crippling loads (C a P^ c 'l? + 2P^ c, f). 
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(psi) 



(lb) 

BLADE 

0.326 

7.60 x 10 6 

0.00309 

23,457 

73,000 

3.11 

1.19 

29,913 

9,100 

1.5 

13,650 

FLANGE 

0.062 

7.47 x 10 6 

0.00309 

23,082 

69,000 

2.99 

1.18 

27,237 

1 .689 

1.0 

1,689 
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Effective— Width from Compressive St r ess Distribution in a 

Buckled Flat Plate 


In order to calculate the crippling strength of the panel, the 
skin, which buckles first, must also be considered. This 
reguires an ef fective- width concept which was originally 
developed for metal structures by T. von Karman (Reference 9). 

In this method, a uniform compressive stress (a '^)# at the same 

average strain as the stiffener at crippling, acts on a width of 

plate w® s directly adjacent to the supported edges. The value of 

w| s is adjusted so the (a Cf ®|) * (w| s ) * (t? k ) is equal to the total 

load carried by the skin on one side of the stiffener. Thus, for 
a skin having the postbuckled distribution shown, the effective- 
width can be found using von Karman 1 s equation. The value of 

(a c,e ?) depends upon the magnitude of the applied design load or, 

in the case of analyzing a tested panel, the failure load. 


es 

W j 



Effective-Width Equation (von Karman): 


w 


es 


= (b?/4)[i + (F;;po7)i 


cr,sk,tt>E c,es 
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The detail buckling and crippling analysis is directed at this 
cross section. The calculations focus on only the middle 
stiffener which carries one-third of the total load up to 
crippling. The predicted crippling load of the middle stiffener 
and skin is equal to the summation of the stiffener element 
pseudo-crippling loads and the effective-width skin load. 


Pj C ;f S - [ CatP^, + 2(P + P C '| S J] 
P c '? s = (a c '| s )*(w| s )*(t? k ) , where wf s = (w? s 

/ x f -L 1 1 1 V 1 

Pertinent dimensions and effective-widths are shown. 
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Load/Strain Curves Across Panel 


Uniform strain was found in the central panels up to skin 
buckling as shown by gauges IN through 4N. Although one of the 
outer panel gauges (IN) is displaced from the others, it has the 
same slope. These gauges indicate that the applied compression 
load was uniform. 


P (10 4 LBS.) 
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Back-To-Back Load/Strain Plots in Inner Skin Panel 


The postbuckling behavior of the inner panels, based upon gauges 
3N & 3F, was moderately nonlinear. This plot indicates that 
buckling occurred between 20,000 Lbs and 25,000 Lbs. 


P (10 4 LBS.) 















Back-to-Back Load/Strain Plots at Tip of Middle Blade 


The lateral buckling of the middle stiffener's blade is indicated 
by the plot of back-to-back gauges 9N and 9F. Initial buckling 
appears to occur at a panel load of about 32,000 lbs, where the 
postbuckling behavior is slight up to a load of about 42,000 lbs. 
Beyond this load level, significant buckling deformation begins. 
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Back-To-Back Load/Strain Plots on Flange of Middle Sti ff P n P r 


The behavior of one of the stiffener flanges is shown by the 
back— to— back gauges ION and 10F. The load— strain plots are 
nearly linear up to a panel load of approximately 40,000 lbs. 
Buckling becomes quite evident at a panel load of about 45,000 

lbs, which is slightly greater than that previously shown for the 
blade . 


P (10 4 LBS.) 
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The compressive load on the stiffened panel is 4b,uuu ids. 
this load, strain gauges indicate that both the skin and blades 
have buckled. Note that the split steel tubes have been mounted 
on the outer unloaded edges. 
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Postbucklinq Behavior of Blade-Stiffened Panel at- 5S r finn T.hs 


At 55,600 lbs., the buckling of the skin, and particularly the 
blades, has become quite severe. However, out of plane 
deformation will become much greater before crippling occurs. 
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ippling load of the stiffened panel was 67,750 pounds. T] 
specimen is shown after being removed from the test rig. 
he severe crimping of the skin and the extensive 
nation of the left blade. The postbuckling forces of the 
skin panels also severely bent the steel split tubes. 
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Stability Analysis Boundary Conditions 


For completeness, buckling and crippling predictions were 
obtained for two boundary condition configurations. 

Configuration A represents simply supported boundary conditions. 
Configuration B represents boundary restraint between simply 
supported and fully fixed. The results are presented in the 
following table. 


Configuration A 

Configuration B 

Skin 

Both Edges Simply Supported 

Skin 

Both Edges Fixed 

Blade 

One Edge Simply Supported 
One Edge Free 

Blade 

One Edge Simply Supported/Fixed 
One Edge Free 

Flange 

One Edge Simply Supported 
One Edge Free 

1 

Flange 

One Edge Simply Supported 
One Edge Free 



Blade ► 

^Flange 


1 

| 


Skin 
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A Comparison of Test and Analysis Results 


This table shows the comparisons between test and analysis for 
the configurations just defined. The analysis is conservative 
for both configurations. However, Configuration B provides much 
closer agreement between test and analysis for both buckling and 
crippling. The predicted buckling strength is about 12 -s ^ 
conservative and the predicted crippling load is about 17 -s 
conservative . 

In conclusion, a test to failure of a blade— stiffened 
carbon/epoxy stiffened panel has been presented. Axial strain 
gauges were employed to verify uniformity of axial strain prior 
to any local buckling. In addition, back-to-back axial strain 
gauges were used for detection of initial buckling and 
postbuckling behavior of the skins, blades, and flanges. The 
stiffened panel behaved as designed. The skins buckled first, 
the blades second, and the flanges last. In the analysis, it was 
assumed that crippling of a blade occurred first, where initial 
failure would be at the supported edge, the location of maximum 
compressive stress. A videotape of the test was made, and it. 
appeared that failure did indeed start at one of the blade/ skin 
intersections . 



Test 

Result 

Analysis 

Configuration A 

Config 

luration B 

Stability 

Mode 

Load 

(Lb.) 

Load 

(Lb.) 

Difference 

(%) 

Load 

(Lb.) 

Difference 

(%) 


20,200 

12 

Buckling 

=23,000 

12,600 

45 

Crippling 

67,750 

41,000 

39 

56,300 
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Symbols and Abbreviations 


b 

b/t 


b i 




b? k 

1 


Ca 


D. . 


E 

F 

F 


13 

c , u 
Th, i 

cc 

cc,u 

* ,i 


F 


cr 



per , sk 
F cl,i 


F cr,u 

* / i 


p cr,u,0E 
cl , i , ss 


area of element "i" , type "u" 
blade 

width over thickness ratio 
Height of blade element "i" 

width of flange element "i" 

width of skin element "i" 

end-fixity coefficient of column: approximately equal to 
3 . 6 for potted end columns in a test machine 

correction factor for edge support of blade and 
stiffener 

flexural/twisting stiffness terms of laminated plate 
in— plane compression modulus of element "i" , type "u" 
crippling stress (psi) 

expected crippling stress of element "i", type "u" 

buckling stress (psi) 
classical buckling stress (psi) 

classical buckling stress of skin element "i" 

expected buckling stress of element "i", type "u" 

classical buckling strength of a long plate type "u" , 
element "i" with simply supported unloaded edges 
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gynihinl s and Abbrevi ations 


F cr f u,0E c i ass i ca i buckling strength of a long plate type "u 
c!,!, x element w ith fixed unloaded edges 


F cr,u,lE classica i buckling strength of a long plate type "u 1 , 
Cl,1 ' SS element "i" with one simply suppported unloaded edge 
pCU ultimate compression strength (psi) 

F cu,u ultimate compression strength of element "i' , type u 

n , i 

and data type "n" 

i element number (i.e., each stiffener has 1 blade, 2 

flanges and 2 skins) 

L column length 

L i effective column length (L f = L * Jc) 

n data type (i.e., empirical "Em", classical "cl", or 

theoretical) 

P C 'J S compression load in effective skin element "i" 

p cc crippling load (lbs) 

p cc ,b expected crippling load of blade element "i" 

* ,1 

P cc ' f expected crippling load of flange element "i" 

* ,i 

cc, ses expected crippling load of stiffener/effective skin 

*, i 

element " 1 " 


p cr buckling load (lbs) 

t thickness (inch) 



u 


thickness of skin element "i" 


element type (i.e., blade b , 
"st", skin "sk", effective skin 
"sts", stiffener/effective skin 


flange 
"es" , 
"ses" 


"f", stringer- 
stiffener/ skin 
panel "p") 
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Symbols and Abbreviations 


w<r s 

1 


w? s 

1 


.cr , st 


r c,es 

, i 


effective width of skin element "i" 

effective width of skin element "i" excluding the width 
of the adjacent stiffener flange 

expected buckling strain of stiffener element "i" 
compression stress of effective skin element "i" 
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REPEATED BUCKLING OF COMPOSITE SHEAR PANELS 


SUMMARY 

Failures in service of aerospace structures and research at the Technion 
Aircraft Structures Laboratory have revealed that repeatedly buckled stiffened 
shear panels might be susceptible to premature fatigue failures. Extensive 
experimental and analytical studies have been performed at Technion on repeated 
buckling, far in excess of initial buckling, for both metal and composite shear 

panels with focus on the influence of the surrounding structure (see for example 
Refs. 1 and 2). 


The core of the experimental investigation consisted of repeated buckling 
postbuckl ing tests on "Wagner beams" in a three-point loading system under 
realistic test conditions. The effects of varying sizes of stiffeners, of the 


magnitude of initial buckling loads, of the panel aspect ratio and of the cyclic 
shearing force, V^, were studied. The cyclic to critical shear buckling ratios, 

(V cyc /V cr ) Were on the hi ® h side > as needed for efficient panel design, yet all 
within possible flight envelopes. The experiments were supplemented by analytical 
and numerical analyses. 


For the metal shear panels the test and numerical results were synthesized in 
Ref. 2 into prediction formulas, which relate the life of the metal shear panels to 
two cyclic load parameters: the (working/buckling) load ratio, (V /V ) , and the 
(ultimate/working) load ratio, (V^/V^l. which reflects the work^g'Ld level 
in the flight envelope, and one geometrical parameter: the (plate/stiffener) 
stiffness ratio, (b 3 t/I f ). It was also found there that the level of shear load, 


at which local yielding first takes place, V dominate 

y 

and hence the life predictions could be expressed in a 


s the endurance of the panel, 
simpler manner, in terms of 


a single load ratio (V /V ) 

eye y 

The composite shear panels studied were hybrid beams with Graphite/Epoxy webs 
bonded to aluminum alloy frames (see Fig. 1). The test results (see Refs. 3 and 4) 
demonstrated that composite panels were less fatigue sensitive than comparable 
metal ones, and that repeated buckling, even when causing extensive damage, did not 
reduce the residual strength by more than 20 percent. All the composite panels 
sustained the specified fatigue life of 250,000 cycles. The extent of damage 

pended on the working load level but no matter how pronounced it was it did 

not affect the fatigue life and did not result in immediate catastrophic failure 
(see Fig. 2 for damages in a typical test). 


*See reference list* 
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The effect of local unstiffened holes on the durability of repeatedly buckled 
shear panels was studied for one series of the metal panels (see Ref. 5). Tests on 
2024 T3 aluminum panels with relatively small unstiffened holes in the center of 
the panels demonstrated premature fatigue failure, compared to panels without 
holes. Even very small holes (of 3 mm diameter and less) caused very significant 
reductions in fatigue life, already at a relatively low load level, for which no 
fatigue failure is predicted in the case of similar unperforated panels. The holes 
caused a shift in the mode of the fatigue failure, initiating now instead of in the 
corners of the shear web in its center (see Fig. 3). Holes with initially 
introduced cracks were compared with smooth ones, the former exhibiting more 
pronounced life degradation, especially for the smaller holes. 

Preliminary tests on two graphite epoxy shear panels with small holes in the 
center showed no similar fatigue life degradation and no shift in failure mode (see 
Fig. 4). Further tests on the effect of holes are in progress. 
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Figure 2. Typical failure mode of a graphi te/epoxy panel (HWB 6A), dynamically 

loaded through two complete "fatigue lives" of 250,000 cycles each, and 
then tested statically from (Weller & Singer, 1990). 


41 




r 



3W0-I3A (jrw) 

~ AS, CLOTH 
D * A 

Vs C y t *l200 


5 %, 

^kshabii 

E3U33I 

raisai 


Figure 3. Fatigue failure of 
2024 T3 aluminum alloy shear 
panel, with a central hole of 
4 mm diameter (WBL 6A), after 
137500 cycles at ^^=1200 kg 

The failure occurs as a crack 
perpendicular to the tension 
diagonal initiating at the 
hole (two initial cracks made 
matters worse here). 


Figure 4. Fatigue failure of 

a graphite/epoxy shear panel 

(HWB 13A), with a central 

hole of 4 mm diameter, after 

250,000 cycles at V =1200 kg 

eye & 

The failure occurs as cracks 
along the tension diagonal, 
emanating from the stress 
concentrations at the 
corners, as in unperforated 
shear panels. 
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Coordinates and geometry of a composite laminate 
with a central circular cutout under compressive loading 


The composites studied are fiber- composite laminate plates made of carbon 
fibers and a thermoplastic-matrix material. The elastic properties of the 
lamina are: E - 15.6 x 10 6 (psi) , E 22 = 0.9 x 10 6 (psi) , u 12 = 0.313, G 12 = 

13 V -\ 10 , (psl >’ and G 23 “ °- 31 ^ 10 6 (psi). The plates have a square 
geometry with a length of 12 (in), a cutout diameter of 2 (in) and a constant 
lamina thickness of 0.005 (in). A [0/90/±45] ns layup is considered. Biaxial 

laminates S aPPll6d ln the f ° rm ° f uniform displacements along the edges of the 
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Solution convergence for transverse shear Q x at (-3., -3.) (in) 
in a clamped [0/90/±45 ] 12s P late without cutout 
under biaxial compression (N x /N y = 2 , t/L = 0.04) 


The transverse shear force Q x is the resultant of r xz integrated over the 
laminate thickness. Q x is interpolated at (-3., -3.) (in) from t e va ues a 
the four Gaussian points of the element containing this location (using a 
bilinear interpolation). Three finite-element meshes are considered. 



0.02 


25 elements 
ooooo 81 elements 
b odo o 169 elements 



Effects of cutout and laminate thickness on maximum shear Q 
in buckling and postbuckling response of a clamped [0/90/±45] ns Opiate 

under biaxial compression 


Without cutout, |Q X max | i s located at (±3.3,0.) for t/L - 0.02 and t/L = 

0.04 and also for t/L =0.06 and t/L =0.08 before activation of higher (i.e. 
second and third lowest) modes takes place for these two thickness/length 

y 1 nc I Kmm XT 1 KT a ' O 


ratios (beyond N= 1.7 N Xcr and N x = 1.5 N Xcr , respectively). After 


# * • ■■ i VV 2{Q £ 7 i “** * *** Jr / * Xi J 

activation of higher modes, the location is at (+6 +4 7) for t/T 

t/L =0.08. ’ f 


0.06 and 


With cutout, |Q X max | is located at (±3. 5, ±1.8) for t/L = 0.02 and t/L 
and for t/L =0.06 before activation of higher modes (N v < 1 7 N ) 

Q f °r t/L = 0.08, I Q x max | is located at the hole free edgfat 
(0 38,_0.92) before activation of higher modes. After activation of higher 
modes for t/L =0.06 and t/L =0.08, the location is at (±6. ,±4.7) . 


0.04, 
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Effects of cutout and laminate thickness on maximum shear Q y 
in buckling and postbuckling response of a clamped [0/90/±45] ns plate 

under biaxial compression 


Without cutout, |Q y max | is located at (0.,±6.) for t/L = 0.02 and t/L - 0.04, 
and also for t/L =0.06 and t/L =0.08 before activation of higher modes takes 
place (beyond N x = 1.7 N Xc r and N x “ 1-5 N xcr , respectively) . After 
activation of higher modes, the location is at (0. , ±4.7) for t/L = 0.06 and 
t/L =0.08. 


With cutout, | Q y max | is located at (0.,±6.) for all four thickness/length 
ratios considered. Activation of higher modes for t/L = 0.06 and t/L = 0.08 
does not change the location of |Qy max l • 
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Effect of mesh refinement on buckling and postbuckling solution convergence 
for a clamped plate [0/90/±45 ] 2 4 s without cutout 
under biaxial compression (N x /Ny = 2, t/L * 0.08) 

For this thick laminate, activation of second and third lowest eigenmodes 
takes place beyond N x - 1 . 5 N Xc r? but no change in buckling mode occurs as the 
structure gradually loses its stiffness and becomes unstable. 
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Effects of cutout and laminate thickness of lowest three eigenvalues 
of a clamped [0/90/±45 ] ns plate under biaxial compression (N x /N y = 2 ) 


The eigenvalue parameter (A N Xo L 2 / D 22 ) is defined in such form that the 
lowest eigenvalue would have the same value for all thickness/length ratios if 
transverse shear was not present. This parameter is plotted with respect to 
the thickness/length ratio. 
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Effects of cutout and laminate thickness on buckling and postbuckling 
response of a clamped [0/90/±45] ns plate under biaxial compression 

(N x /N y =2) 


The load parameter (N x L 2 / D 22 ) is defined in such form that buckling would 
occur at the same value for all thickness/length ratios if transverse shear 
was not present. Likewise, the strain parameter U L / t 2 is such that all 
load/end- shortening curves for the cases with cutout and for the cases without 
cu tout are identical prior to buckling, respectively. 



U L / t 2 



Effect of imperfection sensitivity on transverse shear Q x at (-3., -3. 
in a clamped [0/90/±45 ] ^2s P^ ate without cutout 
under biaxial compression (N x /Ny * 2, t/L = 0.04) 


Three imperfection magnitudes (with respect to the laminate thickness) 
considered: 0.1%, 1% and 10%. The imperfections are made of a linear 
combination of the normalized three lowest eigenmodes. The resulting 
imperfection geometry is close to the first eigenmode (buckling mode) . 



(in) 


are 
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Effect of imperfection sensitivity on buckling and postbuckling response 
(with a change in buckling mode) of a clamped [0/90/±45] ^2s P late 
without cutout under uniaxial compression (Ny * 0, t/L = 0.04) 




Effects of boundary conditions and stress-blaxiality ratio 
on maximum transverse shear Q x in a clamped [0/90/±45 ] ^ 2 S l am ^ nat:e 

without cutout (t/L = 0.04) 





Effects of boundary conditions and stress-biaxiality ratio 
on maximum transverse shear Q y in a clamped [0/90/±45] 12 s laminate 

without cutout (t/L = 0.04) 
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ABSTRACT 

A continuous fiber composite is modelled by a two-element composite 
cylinder in order to predict the elastoplas tic response of the compos- 
ite under a monotonically increasing tensile loading parallel to fib- 
ers. The fibers and matrix are assumed to be elas tic-per fectly plastic 
materials obeying Hill's and Tresca's yield criteria, respectively. The 
present paper investigates the composite behavior when the fibers yield 
prior to the matrix. 


INTRODUCTION 

The elas toplas tic response of fibrous composites has been the subject of a 
number of theoretical s tudies [ 1 -4 ] . When a composite is subjected to uniaxial 
tension loading parallel to the fibers, a two-element composite cylinder (Figure 
1-a) has been frequently utilized to model the composite response. The loading 
direction together with the axisymmetric geometry of the representative volume 
element simplify the mathematical difficulties associated with the equilibrium 
equations. By implementing a traction- free boundary condition to the outermost 
surface of the representative volume element, it becomes possible to construct a 
well-posed boundary value problem when the fibers and the matrix are assumed to 
obey Hill's and Tresca's yield criteria, respectively. A closed form analytical 
solution requires further simplifications such as elastic-perfectly plastic con- 
stituents, perfect interfacial bonding, etc. When the composite constituents are 
assumed to obey the modified yield criterion proposed by Hill[l], the hardening 
effect of the matrix can be taken into account without mathematical difficul- 
ties. However, the present study focuses on a composite with non-hardening con- 
stituents . 



Figure 1. Continuous Fiber Composite 


(a) Representative Volume Element 

(b) Cross section 
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Hill proposed a relatively simple yield criterion which assumes the differ- 
ence between the axial stress and the arithmetic mean value of the radial and 
circumferential stresses to be equal to the yield stress[l], When this yield 
condition is implemented to a composite with elastic fibers surrounded by an 
elastic -perfectly plastic matrix, the entire matrix yields simultaneously, 
resulting in a discontinuity in the slope of the effective stress-strain curve. 

Mulhern, Roger, and Spencer [2] proposed a rigorous analytical solution for a 
two-element composite cylinder under cyclic loading. Their study models an elas- 
tic core fiber surrounded by an elastic-perfectly plastic matrix tube obeying 
Tresca's yield criterion. The resulting composite behavior is almost as bilinear 
as Hill's solution. However, the slope of the effective stress - strain curve is 
continuous because the plastically deformed matrix zone propagates from the 
fiber -matrix interface to the outer surface of the matrix. 

Ebert and Gadd[3] studied a similar problem for an elastic-perfectly plastic 
core fiber surrounded by elastic matrix. Ebert, et al.[4] extended this to an 
elastic-perfectly plastic matrix. However, the application of their numerical 
solution is restricted to a composite in which the Poisson's ratios of the fiber 
and the matrix are equivalent. 

The present paper extends the study of Ebert and other authors [3 ,4] to a 
two-element composite cylinder representing a transversely isotropic fiber sur- 
rounded by an isotropic matrix in which the Poisson's ratios of the core fiber 
and the matrix need not be identical. 


MODEL FORMULATION 

Consider a metal matrix composite reinforced by continuous fibers under uni- 
axial tension loading parallel to the fibers. The globally averaged stress state 
ot the representative volume element is assumed to be one dimensional. The elas- 
tic-plastic response of the bar can be analytically predicted by solving an 
equivalent boundary value problem of a single core fiber which is perfectly 
bonded to the surrounding matrix tube. The volume element representing the 
equivalent boundary value problem is illustrated in Figure 1. The uniaxial ten- 
sion loading in the fiber direction produces a three dimensional stress state in 
oth the core fiber and the surrounding matrix. When the tension loading 
increases monotonically, either the fiber or matrix yields at a certain magni- 
tude of the applied tension loading. Further increment of the tension loading 
induces the yielding everywhere in the composite. The possible yield sequences 
or the composite constituents of the representative volume element can be cate- 
gorized into three cases as shown in Figure 2. The present study provides ana- 
lytical solutions to the first case under monotonically increasing tension load- 

1 ntr ° 



Figure 2. Possible Composite Yield Sequences 
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The entire mathematical formulation of the present study is based on the 
following key assumptions: 

1. The core is assumed to be a transversely isotropic fiber surrounded by an 
isotropic matrix. 

2. Both constituents are assumed to be elastic -perfectly plastic. 

3. The interfacial bonding between the core fiber and the matrix is assumed to 
be perfect throughout deformation. 

4. The core fiber is assumed to obey Hill's yield criterion. 

5. The matrix tube is assumed to obey Tresca's yield criterion. 

6. The axial strain is assumed to be spatially homogeneous. 

Since the geometry of the representative volume element is axisymmetric and 
the loading direction is parallel to the core fiber, the only non- trivial 
equilibrium equation is 


d CJ y- (J ~ (7 Q 

+ = 0 

dr r 


( 1 ) 


The constitutive equations for the transversely isotropic core fiber and the 
isotropic matrix are 
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respectively . 

Since the plastic strain is incompressible, 


(2-a) 


(2-b) 


P P P 
e r + e 9 + e z 


0 


(3) 


It can be shown that Hill's yield criterion becomes the following yield condi- 
tion for the transversely isotropic core fiber under a transversely isotropic 
loading : 




(4-a) 


This mathematical expression is identical to Tresca's yield criterion. 


The surrounding matrix will yield according to one of the following conditions. 


. m m . 

1 °z ' °r 1 

= Y m 

(4-b) 

. m m , 

1 a z - °e 1 

- Y m 

(4-c) 

m m 

1 °z ' °r 1 

■ m m 

" 1 CT z - °B \ " Y m 

(4-d) 

external boundary conditions and interior compatibility are 
o v =* 0 at r = b 

(5-a) 

a r = unique 

at r = a or r = a, c 

(5-b) 

u r =* 0 at r 

= 0 

(5-c) 

u r = unique 

at r = a or r = a , c 

(5-d) 


When both fiber and matrix responses are elastic, the stress state and the 
displacement field in the representative volume element are determined by match- 
ing the radial stress and displacement at the fiber-matrix interface. The prob- 
lem solving procedure for this elastic deformation is simple and straightforward 
as discussed below. b 

Since the stress state in the fiber is always transversely isotropic, 
f f 


f 

CT z = - 2 "LT P + E L e z 


(6-b) 


p ps an unknown constant to be determined, 
the fiber is given by 


The radial displacement in 


u r = C^r 

The strain components become 


(7-a) 


f f 
e r = e 6 = C 1 

From (2-a) , (6) , and (7-b) , 

E 1 = ‘[(l-i'XT ' ^LT^TL.) p / E T + v LT e zl r 
Within the matrix, 



(7-b) 


(8-c) 


( 9 - a) 
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b 2 -a 2 


(9-b) 


m 2 lV m a 
°z = — 3 — rr~ + E m e z 


(9-c) 


u r = & 2 V + ^3/ r 
The strain components become 


( 10 -a) 


■ r u r , r ^2 ' ^3/1 


( 10 -b) 


m o 

€j9 = u r /r = C 2 + C 3 /r^ 

From (5-a) and (5-b) , C 2 and C 3 are determined as functions of P. 
Pa 2 (l+i^ m ) 

c 2 “ 7 , 7 ‘ ^m £ z 

E m ( b2 -a 2 ) 

Pa 2 b 2 (l+i/ m ) 


E m( b ' a ' 

Then the radial displacement in the matrix becomes 


Pa 2 (l+iO 


( 10 -c) 


(H-a) 


( 11 -b) 


u r = — [(l- 2 l / m )r+b z /r] - ( 12 ) 

E m (b 2 -a 2 ) 

The interfacial stress in the radial direction, -P, is determined from 
the displacement compatibility given by (5-d) . 


("m^LT) 


(13-a) 


where 


I-^TT’^LT^Tl] [^ml (l- 2L/ m +b / a ) 


b 2 -a 2 


(13-b) 


It can be shown that the effective axial Young's modulus of the volume element 
is given by 

E c = E L (a/b ) 2 + E m (b 2 -a 2 )/b 2 + 2 (a/b) 2 (»/ m -i/ LT ) 2 /Ap (14) 

Under monotonically increasing axial strain, e z , either the core fiber or the 
surrounding matrix yields first. The yield strains, strengths, and the Poisson's 
ratios of the constituents govern the composite yield sequence. When the applied 
axial strain reaches a certain value, € z y, the core fiber yields first if the 
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longitudinal yield strength of the fiber is much smaller than the yield strength 
of the matrix. The initial yield strain, e z y, a t which the entire fiber yields 
is given by 


e zY = 



1 


E L a 1 + (* / m" l/ LT) (1" 2 i/ Lt) 


(15) 


After the core fiber yields, the surrounding matrix behaves elastically 
until matrix yield at the interface occurs. Since the plastic strain is incom 
pressible and the stress state is transversely isotropic in the core fiber, 
stresses, total strains and plastic strains in the fiber become 


f f f 

a r = a 6 = ' p - a z = ' p + y L (16-a) 


P P P „ 

e r " e 0 ~ - e z/ 2 (16 -b) 

f f f f 

e r = e 9 = u r,r = u r/ r (16-c) 

where -P is the unknown interfacial stress to be determined. 

f f f 

The first strain invariant, e r +e e +e z , together with the uniqueness of the 
radial displacement at the material interface determine the magnitudes of the 
fiber stress components as functions of the applied axial strain. During this 
strain increment where the matrix still deforms elastically, the stress compo- 
nents and the radial displacement in the matrix given by equations (9-a), (9-b) , 

and (9-c) are still valid if the interfacial stress in the radial direction is 
redefined as 
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The radial displacement in the fiber becomes 
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(17 - a) 
(17-b) 
(17-c) 


(18) 


Further increase of the applied axial strain causes yielding of the matrix 
material. The plastically deformed region then propagates toward the outer sur- 
face of the matrix. During this strain increment, the elastically deformed mat- 
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rix and the plastically deformed matrix exist together as shown in Figure 1-b. 
Outside the interface between the plastic matrix and elastic matrix, the stress 
state and deformation are given by the elasticity solution with an unknown^ 
internal pressure, P - *, acting on the interfacial surface between the plastic 
matrix and elastic matrix. Therefore, within the elastically deformed matrix, 
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where P* is determined from Tresca's yield criterion given below. 
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( 21 - a) 


(21-b) 

( 22 ) 


where c is an unknown function of the applied axial strain, e z . 

Within the plastically deformed matrix, Tresca's yield criterion given by 
eqs . (21-a) and (21-b) determines the stress state and displacement field. If 

the Poisson's ratio of the matrix is smaller than that of the fiber, the fiber- 
matrix interfacial stress in the radial direction is always tensile. On the 
other hand, when the Poisson's ratio of the matrix is greater than that of the 
fiber, the fiber-matrix interfacial stress may be either compressive or tensile. 
After the core fiber deforms plastically, the apparent Poisson's ratio of the 
fiber approaches 0.5 as the applied axial strain increases. The interfacial 
stress in the radial direction may be changed from compressive to tensile before 
the initiation of matrix yield. If the matrix yield strain is far greater than 
that of the core fiber, the interfacial stress in the radial direction at the 
onset of matrix yield becomes tensile. Then, from eq . (3) and (4-c) , 
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From eq. (24) and the strain-displacement relationships, the radial and circum- 
ferential stresses are expressed in terms of the radial displacement and its 
gradient with respect to r. 
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The equilibrium equation thus becomes 
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The solution to the above differential equation is given by 
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k = [2(l-„ m )]-V2 

Since the radial stress and displacement must be single-valued at 

r=a, the unknown constants are determined as 

-( 1 + ^m)(E m e z -Y m )b 2 c-( k+ D 

C 2 = ; 

kE m( 1 ” 2l/ m +b ^/ c2 ) 

( 1+ ^m)( E m^2-Y m )b 2 c( k -D 

c 3 = 

kE mO- 2 i/ m +b 2 /c 2 ) 


r=c and 


(28 - a) 


(28-b) 


can also be determined 


The radius of the matrix yield front, 

the applied axial strain, € z , and material properties. However, it is^more^' 
convenient to express the applied axial strain as a function of the radius of 
the n.atrix yield front by satisfying the uniqueness of the radial displacement 

flee of th er '”t t - 1X Dnt11 ^ yUld fr ° nt reaohes the ontermost sur- 

face of the matrix tube, the axial strain is given as 


€ ry 


where 


Y„ 


(l-2u m ) — - (1-2 ^L^)- 




l-2v m +24>2(l+i' m ) - 


»?2 ( 1 ‘ 2i' m ) 


(29 - a) 
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<h = 1 - 


(29-b) 


(b/c) 2 [ (1-^+^/k) (a/c) < k - 1 ) + (l-i/ m - l / m /k) (a/c) - < k+1 > ] 

1 - 2i/ m + (b/c) 2 

(b/c) 2 [(a/c)( k "l)-(a/c)' ( k+ -*-) ]/k 

- 1 - (29-c) 

1 - 2i. m + (b/c) 2 

After the entire matrix yields, it can be shown that eqns . (28-a) and (28 -b) 

should be corrected for further axial strain increment as 


l-(l-2r; 2 /a 2 )(a/b)-( k+1 ) 3 r, 2 


l-2i 


C 2 = 


m 


( E m £ z- Y m) a 3 


aia( k_ n [ (l-2f/ 2 / Q 2 ) (a/b) " 2k - (l-2r/ 2 /a 1 )] 

l-(l-2ry 2 /a 1 )(a/b)( k ' 1 ) 3 i, 2 ] 

E„ 


(30-a) 


l-2i 


Cr, - 


m 


"mJ 


( E m £ z' Y m) " a 3 


Q 2 a* ( k+1 ) [ (l-2r/ 2 /Q 1 ) (a/b) 2k - (l-2t, 2 /a 2 ) 


(30-b) 


E m f u m + ( -*■ ' ) k 1 

1 (l+^ m )(l-2^ m ) 

Em^nT 

“ 2 (l+„ m )(l-2»/ m ) 

a 3 = l-2i/ m - (1-2 i/ lt ) 

E m L E L Y m J 

Then the stress components in the matrix become: 


(30-c) 

(30-d) 

(30-e) 


m 

F f _ y 
iJ m t z 1 m 

+ a 1 C 2 r( k ' 1 ) + a 2 C 3 r'( k+1 ) 

r = 

l-^m 

m 

8 = 

F f -Y 

Lj m t z - L m 

E m [( i , m k+l/2)C 2 r( k - 1 )-( l , m k-l/2)C3r-( k+1 )] 

+ 

l- 2 ^m 

(!+.,„) (l-2./ m ) 


m m 

ff z “ a 6 + Y m 


(31-a) 


(31-b) 


(31-c) 


Further increase of the axial strain, as mentioned in ref. [2], may cause 
another type of plastically deformed matrix region in which the radial and cir- 
cumferential stresses are identical. The present paper, however, does not con- 
sider this case because the infinitesimal strain assumption may not be valid for 
further increase of the applied axial strain. 
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RESULTS AND DISCUSSION 


The effective stress-strain curve for a composite cylinder can be predicted 
by calculating the average value of the axial stress, a z , as a function of the 
applied axial strain, e z , and the mechanical properties of the composite 
constituents. As an example, the effective stress - strain curve of the composite 
studied by Ebert, et al.[4] is demonstrated in Figure 3. The mechanical proper- 
ties of the composite constituents appear in Table 1. In Figure 3, the solid 
lines represent the present analysis. Within the straight line segment (OA) , both 
the core fiber and the surrounding matrix tube are within their elastic limits. 
When the applied axial strain reaches e zY , the entire core fiber yields. The 
next line segment (AB) represents the hardening behavior of the composite with 
plastically deformed core surrounded by an elastic tube. When the applied axial 
strain reaches t z -y in the same figure, the surrounding matrix starts yielding 
from the fiber-matrix interface. This strain can be calculated from eq. (29-a) 
by setting c=a . Then the plastically deformed matrix region propagates outward 
until the entire matrix tube yields. This smooth transient region is represented 
by the line segment BC . The applied axial strain, e z2 , where this transient 
phenomenon terminates can be calculated from the same equation by setting c=b . 
The composite response to further axial strain increase then follows the remain- 
ing line segment. Within the transient region and for the higher value of the 
applied axial strain, the matrix tube material is assumed to be nonhardening 
even though the material hardens signif icantly(Figure 3 in ref. 4). The exper- 
imental results of Ebert, et al.[4] are also plotted in the same figure. 


Table 1. Constituent Properties the Composite [4] 


Material 

Ultimate 

Strength 

(Ksi) 

0.05% 

Offset Yield 
(Ksi) 

Elastic 

Modulus 

(Msi) 

Poisson' s 
ratio 

SAE 4140 
(Core) 

93 

54.9 

28.7 

0.29 

Maraging Steel 
(Matrix Tube) 

318 

288 

25.5 

0.29 


In Figure 4, the radial variations of the radial, circumferential and axial 
stresses in the composite cylinder of which the core volume fraction is 0.5 are 
plotted for two distinctive axial strain values, e zl and e z2 . The stresses in 
the core material decrease slightly as the axial strain increases from e z p. At 
the onset of the initiation of the matrix yield, the axial stress in the matrix 
tube is constant. As the applied axial strain increases beyond t z \, the axial 
in the matrix has its maximum value at the free surface. 
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Stress(psi) Average Axial Stress (psi) 




Figure 3. Composite stress - strain curve 




Figure 4. Stresses in the composite 






CONCLUSIONS 


The present study provides an analytical prediction of the elastoplastic 
response of continuous fiber composites with weaker fibers. The incremental form 
proposed by Ebert, et al.[4] must be replaced by the second order ordinary dif- 
ferential equation given by eq. (26). Furthermore, the present analysis can 
handle the mismatch of the Poisson's ratios as well as transversely isotropic 
fibers. The present analysis will be generalized for the same type of compos- 

^ nder c y clic loading for providing a comparison to the study of Mulhern, et 
a 1 . [ 2 ] 
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Introduction 


The phenomenon of peeling and debonding of thin layers is a subject of interest 
to those concerned with adhesives, thin films, and layered materials. In recent years 
much attention has been focused on such problems as a result of increased interest and 
application of advanced composites and thin film coatings. (See for example ref 1 
^extensive Ust of references pertaining to the subject can be found therein )‘ A 
related problem which is of interest for its own sake but also represents a simple 
example of a tangled adhesive strip and of coplanar delamination interaction is the 
problem of a looped adhesive strip. This is the subject of the present s^. 

rat I?^rf der f h > 6 the problem of an elastic strip which possesses an adherend on 
1 f L k S SU ^ faces ‘ If the stri P ^ deformed so that two portions of 

lot r5 aC H br ° U ^ t int ° contact -. a portion of the strip becomes bonded and a 
ronfi ?- me e lgU £ e We sba11 be interested in determining the equilibrium 

edges a^ puUed apart“ S ' rlP the of the strip when its 

The problem shall be approached as a moving interior boundary problem in the 
calculus of variations with the strip modeled as an inextensible elastica and the 
ond strength characterized by its surface energy.* A Griffith type energy 
criterion shall be employed for debonding, and solutions corresponding to the problem 
of interest obtained. The solution obtained will be seen to predict the intestine 

behavr n ° n w " bond P o mt propagation", as well as the more standard peeling type 
behavior. Numerical results demonstrating the phenomena of interest a^e presented as 

! . J be seen to revea l both stable and unstable propagation of the boundaries 

of the bonded portion of the strip, depending upon the loading conditions 



Figure 1 


*Bottega, W.J.: Peeling and Bond Point Propagation in a Self-Adhered Elastica 

To appear in Quart. J. Mech. and Appl . Math. 
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Formulation of the Problem 


Consider a thin elastic strip which possesses an adherend on one of its surfaces, 
and let the strip be closed on itself in a symmetric manner such that there exists a 
region in the plane of symmetry where the strip is bonded to itself (Figure ). 
Additionally, let the edges of the strip be subjected to equal and opposite forces as 
shown. As a result of the inherent symmetry of the problem, only half of the strip 
need be considered in the ensuing analysis (Figure 2) . The strip thus consists of a 
lifted segment, a bonded segment, and a looped segment. In what follows, all leng 
scales have been normalized with respect to the half length, L, of the entire strip. 

We shall identify each point on the centerline of the strip by its normalized arc 
length s, measured from the edge at which the external force is applied. In so doing, 
the half strip will be divided into four regions; corresponding to the lifted segment 
defined on 0 < s < a, a < s < b corresponding to the bonded segment , with the looped 
segment of the strip divided into two regions, defined by b < s < s* preceding the 
associated inflection point, s*, and s* < s < 1 following the inflection point. We 
shall be interested in assessing the behavior of the above system as a function of the 
magnitude of the applied load or the corresponding edge displacement. 

Let us first define the right handed cartesian coordinate system (x,y) as shown 
in the figure. In addition, let us define the angle 0(s) which measures the angle that 
the tangent of the strip at point "s" makes with the x-axis as s increases (see Figure 
2). One then may easily find the relations 

x(s 0 )-x(s ) = J Sz cos 8 ds and y(s 2 ) -y^) = J** 2 sin 8 ds . (i-a,b) 

Z J. s ^ i 



x 


P 




Figure 2 
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Energy Functional 


rhp A 1 Sha ^ ap ? r ° ach the P roblem as a “°ving interior boundary value problem in 
e calculus of variations, we next define the energy functional (1 given in Table 1 
where U«» corresponds to the normalized strain energy of the segment of the strip 
defined on the domain D i , and is seen to be comprised solely of bending energy as 
the elastica is assumed to be inextensible . The strain energy of the (perfectly) 
bonded segment of the strip, i. e ., the portion of the strip on [a,b] , is thus seen 

tLZTu J r? , Th ^ functional W corresponds to the normaii^d work done by 

the applied load In that expression, P - PL 2 /D corresponds to the normalized 
counterpart of the magnitude of the applied load, P, and D represents the bending 
stiffness of the strip. The functional T corresponds to the "delamination energy" 
bond 6 1 " bY re P r esents the normalized counterpart of the surface energy of the 
bond, y, while a 0 and b 0 correspond to the initial values of a and b respectively 
Finally, we introduce the constraint functional A, with Lagrange multiplier 7 A 
which constrains the deflections of the segments of the strip on D 3 8 and D 4 to be ’ 

rbT"r<l“ “ Ihe fU " Ctl0nS XJ<S * > be^ expressed & of 

A “ * / b cos 6 ds + A S s* cos 0 ds (2e ' ) 

We note that the inclusion of A is equivalent to treating the segments on D and D 
as separate structures and including the work done by the internal force A at s - s* 
may be seen that the line of action of this force must be parallel to the x-axis 
as a result of the support condition at s - 1 (see Figure 2) 


Table 1 


n = r U^-W+r-A ( 1 ) 

where: 

U (i) =J |(d,/ds) 2 ds (2a) 

D. 


0 

III 

(2b) 

a 

W = -P x (o) = P / cos e ds 

(2c) 

o 

r = 2 7 (a-a o ) - 2 7 (b-b Q ) 

(2d) 

« 

5/3 

1 

* 

VI 

,x, 

1! 

< 

(2e) 
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The Governing Equations 

The governing differential equations, boundary conditions, matching conditions and 
transversal ity conditions for the problem of interest are found by invoking the 
principle of stationary potential energy as shown in Table 2 below. In equation ( ) , 
the parameter 6 represents the variational operator. 

The transversality conditions (7a,b,c) result from the moving boundaries during 
peeling or bonding of the strip and thus are associated with equilibrium configurations 
of the system during these processes. 

The intermediate boundaries a, b and s*, as well as the inflection point angles 
a and a* are found as part of the solution to the problem. 

Table 2 


Principle of Stationary P.E.: = 0 

(1) into (3): 

i do. . 

2 ( dP + p i (cos a r cos 9 i> = 0 * (i = 1>3,4) 

e 0 = w/2 


(3) 


(4) 


where: $ ^ = $ 9 * - 8 , « j- «= s(0) > ^ j ^ > 


(5) 


*2,3,4" * ’ a 3,4 

a* = «(s*) 

’ P 3,4 A 


with B.C.s and M.Cs: 


and T.C.s: 


^(a) = « 3 (b) = n/ 2 , 

(6a, b) 

(7a) 

• 4 0) = 0 

(6c) 



« 3 (s*) = « 4 (s*) = a * 

(6d) 

o(b*)-^ 2 i ,. b ^ 

(7b) 

jf cose ds = -jLcos e ds 

u S 

(6e) 

ds, ds. 

— d i - — 4 i =n 

ds s=s* ds s=s* 

(7c) 


71 


Criter ia for Propagation of the "Bond Zone" Boundaries 

The conditions (7a) and (7b) establish the bond zone boundaries during their 
propagation and state that the values of a and b corresponding to equilibrium config- 
urations of the system during propagation of each interior boundary are those for 
which the bending energy densities at the point s = a' and s -= b + are just balanced 
by the energy of the bond. In this context the quantities G(a') and G{b + ) may be 
identified as the "energy release rates" at the bond zone boundaries. The above 
suggests the criteria for propagation of the boundaries of the bonded region of the 
elastica, as listed in Tables 3a and 3b. 

Peeling 

•ii If ’ ^° r SOI ? e lnl ^ i f 1 a “ a o< ec l n - (8a) is satisfied, no peeling will occur and a 
will remain at its initial value a 0 with the lifted segment bending away from the 
plane of symmetry. If, for some initial a = a OI eqn. (8b) is satisfied, the lifted 
segment of the strip peels away from the plane of symmetry such that the value of a 
increases until the corresponding equality (7a) is satisfied. 

Following the above reasoning, we conclude that for the loop to maintain its 
initial configuration, conditions must be such that eqn. (9a) is satisfied. If 
eqn. (9b) holds, the loop would open as a result of excess bending energy at its 
edge, with b taking on smaller values until the energy of deformation is iust 
balanced by the energy of the bond. 


Table 3a 


Criteria for Propagation 
of Bond Zone Boundaries: 


Peeling: 

If G{a‘„} < 2, 

No Peeling 

(8a) 

If G{a „} > 2-i 

Peeling until a 
satisfies equality (7a) 

(8b) 

Similarly, 

If G{b„*} < 2-i 

Loop maintains initial 
configuration 

(9a) 


If G{b 0 + } > 2y Loop grows (b decreases 

until b satisfies equality (7b) (9b) 
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Bond Point Propagation 


If conditions are such that the bending energy G{a } is sufficiently large while 
Gib + ) is sufficiently small, the bond zone boundary "a" will increase while b remains 
at its initial value until a - b'. At this stage, as the resultant bending moment 
m, acting at the bond point satisfies eqn (10) (Table-3b) , and hence acts in a 
clockwise sense, tending to rotate the strip in this sense, while simultaneous y 
there exists a sufficient surplus of bond energy to counter the bending energy of the 
loop and induce bonding at the loop edge, the strip behaves locally as if rolling 
over its counterpart and the "bond point" s = a - b will propagate such that the loop 
closes and shrinks until the corresponding equality (7b) is satisfied. At this point 
the surplus bond energy is depleted and bonding at s = b can no longer occur. 
Equivalently, the growth condition (9b) as well as (8b) will become satisfied and 
conditions will then be such that propagation can occur in both directions simul- 
taneously. Under such conditions, the loop will expire and the strip will separate 

completely . 


Table 3b 

"ROND POINT PROPAGATION ”: 


If (9a) and (8b) satisfied, a - b Q . 


Resultant Bending Moment: 


m ba = [[d./ds]] s = b = [d»/ds] s=b . - [d»/ds] s=b - < 0 (10) 

(thus, m, is clockwise) 

Thus, "bond point" propagates until equality (7b) satisfied 
and loop expires. 
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Basic Integrals 


u^%?rz™T si t * ,<n > ■ 

iS: e : i °" B poin ! ™e es r d ™: £ < ^i ) iSdr“ i ^ lB 

^^pps -x: “*• 

defined'by complete elliptic integrals of the first kind, respectively, 

4> d^ 


F (q i( $) 


f 


/, „ 2 . 2 1 
yi-q . Sin 


and 


W 




(ii-a,b) 


We shall first consider equilibrium 
(^eD 1 ) and of the looped segment (0eD 3 +Di) 
interaction. 


configurations of the lifted segment 
separately, and then examine their 


Table 4 
BASIC INTEGRALS 
Transformation: 

sin(e. /2) = q. sin 0. (li a ) 

qj = sin (otj/2) (i = 1,3,4) (lib) 

(11) into (4), solving for ds and integrating „ 


segment arc lengths t. : 


*i = a 3 [F^i) - F(q r « 1 )/yP~ 

(12a) 

r 3 = S * - h = [F k (q*) - F(q*,t^]/A 

(12b) 

*4= I** - F k (q*)/7I 

(12c) 

where q 3 =q 4 =q* 

(13a) 

*j = sin ’ 1 {l/(q j^)} 

(13b) 


F(q,*) - Elliptic Integral of 1st Kind 
F k (qj) - Complete E.I. of 1st Kind 
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Behavior of the Lifted Region 


The deflection of the edge of the strip at which the load is applied is found 
by solving equation ( 4 a) for d^/ds and substituting the resulting expression into 
equation (i-a), with appropriate limits of integration. Then, upon incorporating 
the transformation (13), we obtain the load edge deflection, A OI as given in Table 5, 
where 

E( qi ,$) = yi-q i 2 sin 2 ^ i d^. and E^q.^ " E(q r */ 2 ) , (iii-a,b) 

correspond to elliptic integrals of the second kind, and complete elliptic integrals 
of the second kind, respectively. The explicit form of the transversal ity condition 
at the "trailing edge” of the bonded region (i.e., at s - a) may be found by solving 
equation (4a) for [d0 x /ds] s _ a and substituting the resulting expression into (7a). 

We then have the condition which (implicitly) defines the location of the trailing 
edge of the "bond zone", during peeling, given by eqn. (15). 

Substitution of equation (15) into equation (12), with i = 1, and equation (14) 
gives explicit relations for the magnitude of the applied load as a function of "a" 
and the normalized load point deflection as a function of a respectively , with the 
load point rotation a a parameter. Specific results corresponding to selected values 
of 7 will be presented in a later section. 

Table 5 

BEHAVIOR OF LIFTED REGION 


LOAD POINT DEFLECTION: 




-[ 2 E k ( qi )-F k ( qi )]} (14) 


where E(q,*) - Elliptical Integral of 2nd Kind 


E k (q j) - Complete E.I. of 2nd Kind 


T.C. @ s = a" (7a) becomes 

P = - 27 /cos a (jt/2 < a < ») (15) 
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Equilibrium C onfigurations of the Loop 

Th f angle which measures the rotation of the tangent of the loop at its 
inflection point s - s* is found by imposing the condition (6e) . Thus, solving 

left^nc^rieht haV^d ’ substitutin 8 the resulting expressions into the 

III*? t ™ l S ° f , ( , 6e) respectively, noting (6d) and incorporating the 

transformation (H) for i - 3,4 we obtain the condition given by (16) below where 

H6? fhat r r h a “ 3 V\ and q * " q 3 “ <U- Ic ma y be seen from equation 
(16) that a* is independent of the size of the loop, of the material and geometric 

properties of the strip and bond, and of the magnitude of the applied load and 
thus is a "characteristic angle" of the problem. Equation (16) may be solved 
numerically to yield the value of a* as given below. 

t ,-^ he i t0ta u arc len & th of the loop, i, is simply comprised of the sum 

°! f , f” gth f ° f * tS const rtuent segments. Thus, adding equations (12b) and 
(12c) yields the relation for i, given by (18). The relative portions of the loop 
corresponding to its constituent segments are then found by dividing eqns (12b') 
and (12c) by (18). We thus have 7 8 eqns 


and 


V* - t F k<q*) - F(q*,$*)]/[2F k (q*) - F(q*,**)] 
ZJZ - F k (q*)/[2F k (q*) - F(q*,$*)] , 


(19a') 
(19b' ) 


which are seen to correspond to "characteristic length ratios" 
above ratios may be evaluated, using the computed value of a* 
given at the bottom of Table 6. 


of the problem. The 
to yield the values 


Table 6 

EQ UILIBRIUM CONFIGURATIONS OF THF. I OOP 
Imposing (6e): 

2[2E ]< (q*) - F^q*)] = 2E(q*, #*) - F(q*, #*) (16) 

A. A 

(where #* = #^ = 

Solving (16) yields "characteristic angle" of inflection point 

a* = 117.54" (for any loop size and 

mat’l./geom. props.) (17) 

(12b) + (12c) * loop (half) arc length t. 

1 = *3 + *4 = t 2 F |/ 9 *) ' F (q*. **)]/ 7 i" (18) 

(12b,c)/(18) = "Characteristic Arc Length Ratios": 


76 


= 0.3254 and t^/i = 0.6746 
(for any loop size and mat’l./geom. props.) 


(19a,b) 



Debonding of the Looped Segment 


We may next consider equilibrium configurations of the looped segment of the strip 
as it opens (debonds) by evaluating the explicit form of the transversality condition 
at s » b + in a manner analogous to that done earlier for the lifted segment. Doing so 
we find that during opening of the loop, the condition (7b) takes the form given by 
eq'n.(20) where a * is given by (17). Since, as discussed earlier, a* is a "character- 
istic angle" (i.e. , it maintains a fixed value for any equilibrium configuration of the 
loop,) it is seen from the above expression that the internal force A maintains a 
constant value during debonding of the loop. 

Substitution of (20) into equation (18) yields a critical value of the loop length 
given by equation (21) below, with the inequalities (9a, b) now interpreted in terms of 
the arc length of the loop; e.g. - (22a,b). 

It may be noted that a minimum value of the normalized bond energy is required for 
the elastica to remain adhered to itself. This value 7 - 7 min corresponds to the 
limiting case where the loop traverses the entire strip and is in self-contact only at 
the loading edge s - 0 (i.e., it corresponds to the limiting case when 2 cr - 1). Upon 
employing the result (19) we find the desired value given below. Adherends whose 
normalized bond energies possess magnitudes which are below this value are thus not 
"strong enough" to maintain a self-adhered configuration. 

Table 7 

T.C. @ s = b + (7b) becomes: 


x = -2 7 /cos a* (= constant for given 7) (20) 


(20) into (18): 


* = [2F j^q*) - F(q*, t*)]y(-cosa*)/2 7 (21) 


If t > t 

cr 

No debonding of loop occurs 

(22a) 

If t < t 

cr 

Debonding of loop occurs 
t increases (b decreases) until t = 

(22b) 

t = 1 => 
cr 

7 = 7 • = 2.292 

' min 

(23) 
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Effective Bond Strength and Propagation Behavior 


A plot of 2 CI versus 7 is displayed in Figure 3. It may be observed from the 

igure that the amount of "effective bond strength" gained, as measured by the relative 
decrease in 9 c an f nanfiv — . , J Live 


t O O t Hiuaouicu 

significantly decreases as 7 is increased beyond 200 . 


Peeling and Bond Point Propagation 

The above solution offers the following scenario for a looped adhesive strip with 
given 7 > 7 min existing in an initial configuration such that 2 - 1 -b >2 for 
equivalently b Q < b cr - l-i cr ) , with an initial lift zone size o°f 2 *= °a As P is 
increased the corresponding value of a, = a is increased according to equation ( 12 a), 
with the associated deflection A c varied according to equation (14). This process is 
continued until equation (15) is satisfied at which point peeling begins with the 
trailing edge’ of the "bond zone" s = a propagating (and i x increasing). As the 

Jhe in d P g ° 1S . larg6r than the critical length, the loop edge boundary of 
the bond zone remains at its initial value until a - b c . At this point if conditions 
are such (equation (15) satisfied) that peeling continues, the bond point a = b 
propagates with 2 decreasing (b increasing). During this phase the values of a* 

2^/2, and 2J2 maintain the values given by (17) and (19a, b) respectively. The loop 
thus shrinks in size during this phase, with its geometry evolving through successive 
self-similar shapes (as if the strip were being pulled through a rigid clatp at s - b) 
This process continues until 2 - 2 at which point conditions are such that peeling 

tMs°r U r *1 rh th i S " b and S " b (i e > in both directions) simultaneously. At 
this instant the loop is terminated and the surfaces on each side of the plane of 

symmetry separate. Results corresponding to specific values of 7 are presented 



1000 


Fig. 3. Variation of Critical Loop (Half) Length with Normalized Bond Energy 
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Numerical Results 


Results are presented corresponding to selected values of the normalized bond 
energy. Specifically, we consider the strip/adherend systems whose material and 
geometric properties are characterized by the values 7 — 5, 10, 50, 100, and 1000 
respectively. 

The projections of the associated lift zone/bond point "propagation paths" 
in the A 0 -a, P-a and P-A 0 spaces are calculated and are displayed in Figures 4-7. 

Each curve is terminated at the critical values a — b cr = l-i cr , which are given 
by the values b cr - .3230, .5213, .7859, .8486, and .9522 for the respective values 
of 7 considered. The prepropagation load- deflection behavior of lift zone segments 
corresponding to initial lengths of a 0 = 0.25 and a Q - 0.50 are also displayed in 
Figure 6 . Finally, the variation of the magnitude of the internal force A , as a 
function of the loop length i, is displayed in Figure 8 . 

It may be seen from the figures that propagation of the lift zone boundary or 
bond point occurs in a stable manner for a deflection controlled test, and. in an 
unstable and "catastrophic" manner for a force controlled test. The following example 
illustrates the general behavior of the self- adhered elastica. 



Fig. 4. Lift Zone/Bond Point Propagation Paths (7 - 5, 10, 50, 100, 1000): 
Load Edge Deflection vs. Bond Zone Boundary 
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Let us consider a strlp/adherend combination with 7 - 50 which is initially 
configured such that a_ - 0.25 and b„ - 0.60. It may be observed that a system 

site ra tMs rl sLn y ( h ne th°, T S T Ue a r ° f 7 COnsld « ed "> uld "« “Infin a loop 

size this small (b this large) and that for such a system the loop edge boundary b 

would immediately decrease to its corresponding critical value. Returning to the ’ 

previousiy defined case (7 - 50) , the system initially follows the prepropagation 

path for.. - 0.25 in Figure 6 . as P and 4„ increase from the origi„ P During this 

phase the system simultaneously follows a purely vertical path, corresponding to 

?„ , the 4 -;?, and P ' a =P ac * s (Figures 4 and 5). The system continues to behave 

in this fashion until the propagation path corresponding to 7 - 50 is interceDted 
At this point the lifted segment of the strip has accumulated enough bending energy at 
the bond zone boundary s -a', for the lift zone to propagate. We" note thft as ^ 
*o> *cr (b o < b er ) the bending energy of the looped segment at s - b + is insufficient 
for propagation so b remains at its initial value. Let us first consider the case 
where the load edge deflection, A c , is controlled. For this case, as A is 
incrementally increased, a increases incrementally following the corresponding path 
in Figure 4. The corresponding values of P may be observed 8 from Figs . 5 and 6 to 
decrease accordingly. This process continues, with the strip peeling from its 
symmetric counterpart in a stable manner, until a - b Q . At this point the bending 
whnf that ^ loop at s - b+ is still sufficiently low as to maintain the bond 6 

further ^heV H • 66 6 T 6h t0 Permlt debondin g- As A 0 is increased 

further, the bond point s - a - b then propagates in a stable manner, with the loop 

shrinking through a series of self-similar shapes until i - £ at which point ? 

sufficient bending energy exists on both sides of the bond point and the strip 

separates . For the case where P rather than 4„ Is controlled, the system behaves 



Fig. 5. Lift Zone/Bond Point Propagation Paths (7 - 5 , 10, 50, 100): 
Applied Load vs. Bond Zone Boundary 
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in an analogous manner except that all behavior during the propagation phases occur in 
an unstable manner. Thus, for this case the process of complete separation of the 
elastica is initiated as soon as P reaches a critical value. 

The phenomena described above may be observed by simply closing a piece of 
adhesive tape on itself, thus forming a loop, and then peeling the edges apart. Such 
an "experiment 1 ' would correspond to a deflection controlled test, with the normalized 
bond energy characterizing the tape observed to be at the upper end of the range 
of values considered in the numerical simulations presented herein. 



Fig. 6 . Prepropagation and Lift Zone/Bond Point Propagation Paths (7 5, 10 

50, 100): Applied Load vs. Load Edge Deflection 
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BASIC IDEA 


Due *° ,^ rtmn 1 fype Of loading and geometrical boundary conditions, each beam will 
respond differently depending on its geometrical form of the cross section and its 
material definition. As an example, consider an isotropic rectangular beam under pure 
bending Plane sections perpendicular to the longitudinal axis of the beam will remain 
plane and perpendicular to the deformed axis after deformation. However, due to the 
Poisson effect, particles in the planes will move relative to each other resulting in a form of 
anoclastic deformation. In other words, even in pure bending of an isotropic beam, each 
cross section will deform in the plane. V 

^ 1,16 be f n J al * >ve is replaced by a generally anisotropic material, then the 
cross sections will not only deform in the plane, but also out of plane. Hence, in general 
both in-plane deformations and out-of-plane warping will exist and depend on the 
geometrical form and material definition of the cross sections and also on the loadings 
For the purpose of explanation, an analogy is made. The geometrical forms of the bodies 
of each individuals are unique. Hence, different sizes of clothes are needed. Finding the 
sizes of clothes for individuals is like determining the warping functions in beams. 

A new beam theory using first-order warping functions is introduced. Numerical 
examples will be presented for an isotropic beam with rectangular cross section The 
theory can be extended for composite beams. (Fig. 1 .) 



Figure 1 . Analogy between determining the (first-order) warping functions 
and the proper size of clothes for individuals. 
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CANTILEVER BEAM 


Consider the case of an isotropic rectangular cantilever beam with a tip loading ( P ). For 
the purpose of comparison to the St. Venant elasticity solution, St. Venant boundary 
stresses shall be taken into account. These self-equilibrated boundary stresses are shown in 
the figure below. XYZ is the system coordinate ; L is the length ; H is the height, and B is 
the thickness of the beam. Comparison will be made with respect to the plane stress St. 
Venant elasticity solution. Hence, the comparison will be more valid as the thickness goes 
to zero. (Fig. 2.) 



Figure 2. Cantilever beam with its boundary conditions. 


ASSUMED DISPLACEMENT FIELD 

The u,v,w are displacement components parallel to the x,y,z coordinate axis (refer to figure 
2.). -uo is the transverse displacement of the axis of the beam parallel to the x coordinate 
axis. 9 is the bending rotation of the beam axis and its positive sense is defined as in the 
direction of the positive y axis. M and Q are the bending moment and shear force, 
respectively, q is the distributed load. U and V are the in-plane deformation functions 
parallel to the x and y coordinate axis, respectively. W is the out-of-plane warping function 
parallel to z coordinate axis. E is the Young modulus, and I is the moment of inertia of the 
cross section. 

Strains are computed from the first set of equations. From Hooke's law [ (o) = [c](e) ], 
stresses can be calculated. Using the definition of moment, M can be solved in terms of E, 

I and 9. The final form of the displacement field is as shown in the second set of equations 
The detail is as follows. 

By definition of moment, 

M(z) « J xozdA 

Assuming that (in consistency with using only the first-order warping functions) 

Q(z) = - q( Z ) = o, the moment will be expressible as being M W = -E I <pVz). By equilibrium 

Q(z) = m’(z) =- E i (ffrz). Substituting M and Q into the first set of equations , the final form of 
the displacement field is obtained. 

It is important to note that no assumption is being made except the assumed displacement 
field itself. (Fig. 3.) 


u(x,y,z) = u 0 (z) + M(z) U(x,y) 
w(x,y,z) = - x <p(z) + Q(z) W(x,y) 

where 

U(x,y) = - (x 2 - y 2 ) 

W (x,y) = W (x) = - x3 
6 E I 


Final Model 

u(x,y,z) = u 0 (z) + <p'(z) U(x,y) 

w(x,y,z) = - x <p(z) + tp'(z) W(x) 
where now 

U(x,y)= y (x2-y2) 

WTt \ 2 + 0 

W(x) = — - — x 3 
6 


v(x,y,z) = M(z) V(x,y) 


o 


V(x,y) = ---xy 


v(x,y,z) = (p (z) V(x,y) 


V(x,y) = o x y 


Figure 3 . Proposed model in a case of rectangular cross section. 



LAYOUT OF THE NODAL POINTS 


The finite-element model is developed using a layout of the nodal points as shown below. 
The layout is chosen such that all terms in the strain energy expression are taken into 
account. The minimum order of polynomials that is required for v based upon the strain 
energy expression is three. Hence a four-node layout is required for <p. A five-node layout 
is selected for u o because, from the physical point of view, uo is one order higher. Other 
polynomials could also be selected (Fig. 4.) 




uo 



Displacement in the X-direction 


Bending rotation 


Figure 4. Layout of the nodal points in the finite-element model. 
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FINITE-ELEMENT MODEL 

Using the principle of stationary potential energy a finite-element model is obtained. The 

term (sm } in the equation is due to the fact that, instead of applied concentrated resultant 
forces, 'applied' St. Venant distributed stresses (through the cross section) are to be 
considered. If applied concentrated resultant forces do exist in the reality, then this term 

will vanish. The term { 5m ) is present due to the fact that the distributed forces are applied 
on the upper surface of the beam. These additional terms exist because a beam theory that 
accounts for in-plane deformations and out-of-plane warping is used. Had an Euler- 
Bemoulli beam theory been used (or likewise Timoshenko beam theory), all these terms 
will vanish no matter how the loads are applied. All other terms are the usual terms that 
result when developing a finite-element model based upon an Euler- Bernoulli theory. For 
example 

(f}=fq[N] T dz 

where q is the distributed load and [ N ] are the shape functions. (Fig. 5.) 


5 

Uo(z) = luoj 4>Xz) 

j=l 

4 

<p( z ) = Z<pj e / z ) 

j=l 

Uoj , <Pi are Lagrangian type shape functions. 

Equilibrium equations 

[Kll] 5 x5 [K12] 5x4 1 / { Uo >5x1 | / { P ) 5 x 1 ’ 

[K2l] 4 x 5 [K22] 4 x4 \{tp}4xlj ^ {M } 4 x 1 ) 

{0}5xl \ / {f}sxl 

+ + 

i (6m) 4x1 ) \ (8m) 4x i / 


Figure 5. Finite -element model. 
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TIP DISPLACEMENT OF A CANTILEVER BEAM 

A semi-logarithmic plot of aspect ratio versus nondimensionalized tip displacement is 
shown below where the tip displacement is nondimensionalized by dividing by the Euler- 
Bemoulli solution for a given load and geometry. Aspect ratio is defined as the length 
divided by the height of the beam. The Poisson ratio is taken equal to 0.25. All theories are in 
agreement for slender beams. For this type of loading, the elasticity solution, the proposed 
theory, and the Timoshenko using k equal to 2/3 are in perfect agreement. Using a k factor 
equal to 0.8475 [ 1 ] in the Timoshenko theory results in a stiffer beam (compare to using 

k equal to 2/3). . , 

Solutions were calculated for extremely long slender beams ( IVH - 100 ) to insure tnat 
the current beam element converged to the Euler-Bemoulli solution and did not shear- 
lock”. All calculations were performed assuming that B/L is equal to 1/8000. This 
selection was made to insure that the current model can be directly compared to the plane 
stress elasticity solutions. (Fig. 6.) 



Aspect Ratio | IVH I 


Figure 6. Tip displacements of a tip-loaded cantilever beam. 
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[ ( Z/H ) / X ] 


NORMAL STRESSES AT THE ROOT OF A CANTILEVER BEAM 

Consider a case of the cantilever beam with an aspect ratio equal to three. The abscissa is 
the nondimensionahzed normal stresses with respect to the Euler-Bemoulli normal stress at 
the top of the surface, i.e., . The ordinate is the nondimensionalized X-coordinate where 

the top and bottom surfaces of the beam are defined as -1 and 1, respectively As can be 
seen, all theories are in perfect agreement. (Fig. 7.) 



[ CTz / ( O z )Eu!er-Bemoulli ] 


Figure 7. Normal stress at the root of a tip-loaded cantilever beam. 
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TRANSVERSE SHEAR STRESSES 
AT THE ROOT OF A CANTILEVER BEAM 


Nondimensionalized shear stresses with respect to the elasticity shear stress at x equal to 
zero are made, i.e., nondimensionalized with respect to ta. The ordinate is the 
nondimensionalized X-coordinate. The results from the proposed theory are in perfect 
agreement with the elasticity solutions. If the Timoshenko theory is applied constant shear 
stress distribution is obtained. In fact, their values are equal to the average shear stress, i.e., 
P/A where A is the area of the cross section and P is the applied concentrated load. In t 
case they do not depend on the value of the shear correction factor (k). Hence, the shear 
correction factor will influence the stiffness and shear strains of the beam, but not the shear 
stresses Since the shear stresses are independent of k, then the shear strains must vary 
proportionally to the inverse of k. For an Euler-Bemoulli beam a contradiction exists It 
shear stresses are computed from the shear strains, then their values will vanish On the 
other hand, from the equilibrium point of view, shear stresses cannot be zero. Using the 
principle of equilibrium, shear stresses can be obtained. (Fig. 8.) 



0.0 0.2 


Elasticity 


- Timoshenko (k independent) 


[ Tzx / ( Tzx plasticity ] 


Figure 8. Transverse shear stresses at the root of a cantilever beam. 


TABLE OF COMPARISON FOR THE STRESS COMPONENTS 
OF A CANTILEVER BEAM 


Results from the theory of elasticity for the stresses are available. Due to the computational 
round-off errors, the coefficients e i will exist in the proposed theory. These terms become 
smaller as the value of B (the thickness) goes to zero. If truncation errors could be 
eliminated, the E i terms will go to zero as B goes to zero. This is due to the fact that the 
out-of-plane warping function W(x,y) was taken from the plane stress solution, which is 
then only a function of x, i.e., W = W(x). This is done mainly for the purpose of 
comparison and simplicity. In order for the current comparison to be valid, the thickness 
of the beam should be taken very small ( B --> 0 ). As one can see, apart from the round- 
off errors, the proposed theory is in perfect agreement with the elasticity solution for the 
whole body (plane). (Fig. 9.) 


Oz = ei X z2 + CiJLZ + £2 x3 + e 3 X 

Tzx = £4 Z3 + £5 Z2 +£6 X2 Z + £7 Z + C 2 X 2 + £3 
(Tx CTy — £g X3 Tyi — Txy = 0 
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Elasticity 
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Elasticity 

Cl 
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1.00 

1.00 
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0.00 

0.00 

C3 

1.00 

0.67 

0.00 

£i — > 0 as B (the thickness) goes to 0 


Figure 9. Comparison of the stress components relative to 
the elasticity solutions. 




SIMPLY SUPPORTED BEAM 


Since the model is developed for tip loadings, it may have slight deficiencies for 
analyzing beams with higher-order distributed loads. In this example, the effects of the 
higher-order warping functions will be shown. Although it is possible to extend the model 
incorporating some degree of higher-order warping functions, for beam-type structures 
it may not be necessary. Higher-order warping functions play an important role for beams 
with small aspect ratios (closer to solid structures). But, as the aspect ratio gets smaller, 
the St. Venant solution becomes trivial in the practical sense (in reality). Another way of 
defining beams is as follows. A structure (structural element) can be considered as being a 
beam if the higher-order warping functions play insignificant roles. (Fig. 10.) 


m q 





Figure 10. Simply supported beam with its boundary conditions. 



MID-LENGTH DISPLACEMENT OF A SIMPLY SUPPORTED BEAM 


It is very interesting to note that, although the Timoshenko theory for k equal to 2/3 gives 
the exact results (for the displacements) in the case of a cantilever beam, it gives the most 
flexible structure in the case of simply supported beams with constant distributed load. The 
proposed theory still gives very accurate results for the displacements. As can be seen, the 
effect of higher-order warping functions in the displacements, in this typical case, is very 
insignificant. (Fig. 11.) The effects in the stresses can be seen in the next figure. 



Aspect Ratio [ L/H ] 


Figure 1 1 . Displacements at mid-length of a simply supported beam. 




NORMAL STRESSES AT MID-LENGTH OF 
A SIMPLY SUPPORTED BEAM 


Consider the case of a simply supported beam with a uniformly distributed lo ^ d ^ ere 
the aspect ratio (L/B) is equal to three. Nondimensionalized normal stresses are with 
respectto th^Eule^Bemoulli stress on the bottom surface of the beam.. It is important to 
note that the elasticity solution results in a cubic normal stress distribution as opposed to linear 
as suggested by the Euler-Bemoulli or Timoshenko beam theories. (Fig. 12.) In the next 
figure, a closer look at the lower portion of the cross section is presented. 



Figure 12. Normal stresses at mid-length of a simply supported bean . 
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NORMAL STRESSES AT MID-LENGTH FOR 
THE LOWER PORTION OF THE CROSS SECTION 


rSth Ca H bC SCen ’ the elasticity solution gives higher normal stresses at the top and bottom 
of the beam compared to the Euler-Bemoulli or Timoshenko beam. For this typical 

EltS (V 1 th an as P ect L ratio equal to three), the elasticity normal stiess at the 
f th< V° P 1S 3 Percent higher, and the proposed theory gives 2.7 percent higher 
This effect is due to the presence of distributed load or, in other words, due to the presence 

signilc?nt 0r (Fig W ^) nS functlons ' As the as P ect ratio gets smaller, its effect will be more 



Figure 13. Detail of figure 12 for lower portion of the cross section. 



TRANSVERSE SHEAR STRESSES AT LEFT-END 
OF A SIMPLY SUPPORTED BEAM 


The same nondimensionalization as is found in figure 8 is made. The proposed theory 
results in perfect agreement with the elasticity solution. Again, the Timonshenko beam gives 
a constant shear stress distribution which is equal to R/A where R is the reaction force at the 
left end. (Fig. 14.) 



Figure 14. Shear stresses at the left-end of a simply supported beam. 



TABLE OF COMPARISONS FOR THE STRESS COMPONENTS 
OF A SIMPLY SUPPORTED BEAM 

All underlined terms exist in the plane stress elasticity solution. The terms with the 
coefficients hi and 1\2 are the nonclassical terms and become important only for beams with 
very small aspect ratios. As a matter of fact, these two terms, found in the expression of c z , 
are self-equilibrated in the section planes. Obviously both the Timoshenko and Euler- 
Bemoulli beam theories cannot capture these higher-order terms. The proposed theory is still 
able to capture these two terms. Although it is not accurate inside the body of the beam, it 
gives an accurate result at the top and bottom of the beam, which are, in fact, the most 
important points (for normal stresses) for practical purposes. 

Due to the presence of the distributed load, the stress component a x will not vanish. The 
proposed theory gives a meaningless result in the sense that it does not satisfy the boundary 
conditions at the top and bottom of the beam. This is not surprising because the theory accounts 
for only first-order warping functions. Although it is possible to extend the theory 
incorporating higher-order warping functions, it may not be necessary for practical 
purposes. It is important to note that for this typical beam structure, the most important stress 
components are <J z andT zx . (Fig. 15.) 


Cz = Ci x 2. 2 + ClXZ + hj_x3 + h? x 

Izx = El z3 + e 2 7.2 + Cr X2.Z + £4_Z + £s_X2 + £6 

Ox = C7JC3 + C8_X + £9 Oy = d] X3 Xyz = Txy = 0 
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Figure 15. Comparison of the stress compoments relative to the 
elasticity solutions of a simply supported beam. 




FUTURE STUDIES 


Future studies will include vibration analysis, methods to determine the warping functions 
for a typical beam cross section, geometrical nonlinearity, and an extension for composite 
beams. Basically the proposed theory will be extended for general cross sections and 
material definition (composite beams) such that they can be applied for any special case. 
(Fig. 16.) 


FUTURE STUDIES 
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Figure 16. Future studies. 
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STRUCTURAL MODELS FOR SYSTEM RELIABILITY EVALUATION 

E va l U ati° n of structural system reliability is based on the determination of 
probabilistic response of structures. The study of probabilistic behavior of 
structural systems requires an underlying structural analysis model. The current 

S : m n el r bl lty / nal r iS procedures utilize the structural analysis methods 
methnd lly h 6V a l0 ? ed Wlt \ a dete rministic point of view. The applicability of these 
methods when dealing with a random situation needs examination. It may be that 

methods of structural analysis that are quite suitable for deterministic analysis 
are not as suitable for probabilistic analysis and vice versa. Development of 
structurai system models especially suited for random state variables may be more 

A recent NaM ^ / t0 f t ° im P ortant insights into random structural behavior. 
A recent National Science Foundation Workshop on Structural System Reliability held 
at the University of Colorado, Boulder [1] also emphasized the need for the 
development of structural analysis models from the probabilistic viewpoint. 

The objective of this paper is to demonstrate discrete extremum methods of 
structural analysis as a tool for structural system reliability evaluation 
r“ 3l y ’ l lneaV and multiobj ective linear programming models for analysis of 
rigid plastic frames under proportional and multiparametric loadings, respectively 
are considered Kinematic and static approaches for analysis formV primal -^II 
pair in each of these models and have a polyhedral format Duality relations link 
extreme points and hyperplanes of these polyhedra and lead naturally to dual 
methods for system reliability evaluation. 3 


m DIRECT METHOD 

- Formulation of the global limit state surface 

Computation of probability of random variables 
having an outcome in the safe set 


■ INDIRECT METHOD 

Identification of all failure modes 

Computation of the probability of failure of 
individual modes 

Evaluation of system reliability from modal 
probabilities 
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LIMITATIONS OF CURRENT APPROACHES 


It is known that the structural problems can be analyzed by idealizing these 
mathematically as either continuous or discrete models, leading respectively to the 
solution of differential or algebraic equations/inequalities. The underlying 
solution principles can be based on the solution of simultaneous equations or the 
use of extremum principles. In the sixties, Schmit and coworkers [2] briefly 
explored extremum methods for the deterministic structural analysis problems but 
found them to be not competitive with solution methods for simultaneous equations. 
However, in the use of simultaneous equations procedure for probabilistic 
structural analysis, one has to solve the structural analysis problem repeatedly 
for different realizations of random variables and this is computionally costly. 

Use of extremum principles, on the other hand, elucidates the mathematical 
structure of the problem corresponding to various random realizations of state 
variables. This structure is extremely coherent with a definite pattern about the 
solutions of the problem. An understanding of such patterns leads one to gain 
important insights into response under random variables without analyzing the 
structure for all such combinations. This coupled with use of recent computational 
advances in algorithms [3] and vector processing of information on supercomputers 
are likely to make these methods extremely attractive for use in probabilistic 
analysis. For example, recent research [4] shows extremum methods to be ideally 
suited for structural analysis required in the system reliability assessment of 
structures with rigid plastic material behavior. 


■ STRUCTURAL ANALYSIS FOR RELIABILITY EVALUATION 

Discrete/Continuous models 

Classical methods, elastic/plastic 
Extremum methods 

■ RESPONSE PATTERNS 

Polyhedral response regions 

Other response regions 


■ ADVANTAGES OF EXTREMUM METHODS 
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EXTREMUM METHODS FOR STRUCTURAL ANALYSIS 

Structural analysis requires both the calculation of the distribution of forces 
throughout the structure and the displaced state of the system under the action of 
applied loads. One of the fundamental features of structural analysis is the 
possibility of using either forces or displacements as basic variables, with the 
respective approaches referred to as static and kinematic methods. The algebraic 
relationship of the static and kinematic approaches are the mathematical transpose 
of each other, a feature known as the static-kinematic duality. If the structure 
is statically determinate, the number of equations is same as the number of 
variables, and the forces and displacements can be found easily from the solution 
of system of algebraic equations. For statically indeterminate structures, 
additional conditions reflecting the material constitutive relations must be 
introduced. 


As an example, for redundant frame structures, partial satisfaction of structural 
constraints generates a subspace in R n containing the solutions of interest and the 
final solution can be reached by some optimality criterion. The optimal solution 
gives the result corresponding to the solution by traditional methods. The power 
of the extremum methods, however, is that all the suboptimal solutions may also be 
obtained from the model, and these suboptimal solutions correspond to various 
random realizations of the variables. This set of available solutions has a rich 
underlying mathematical structure and such patterns have recently been studied for 
rigid plastic frames under proportional and multiparametric loading [4]. 
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RIGID PLASTIC FRAMES WITH PROPORTIONAL LOADING 


The problem of limit analysis of frames in which the plastic behavior is activated 
by a single stress resultant (such as flexure) may be formulated in terms of a 
Linear Programming (LP) model [5]. Plastic hinges are assumed to form at a set of 
discrete locations (j = 1, 2, . . . , J) and the plastic moment capacity at the j t ^ 1 

section is denoted by Mpj . Models formulated from dual structural consideration of 
static and kinematic principles have been shown to be a primal-dual pair in the LP 
format. The variables are Mj = moment at section j , Mj — Mj + - Mj ~ where Mj + 
and Mj represent the positive and negative values of moments; 0-? = rotation at 
section j ; t^ = a coefficient indicating the contribution of the k elementary 
mode to collapse, t^ = t^ + - t^ and X + , X = collapse load factors for the 
kinematic and static LP’s respectively. The parameters ared^-j = hinge rotation of 
elementary mechanism k at joint j , = external work associated with elementary 

mechanism k, M = number of equations of equilibrium and Mpj = member capacity at 
section j . 


KINEMATIC LP 


STATIC LP 


X= Min [j M p (0,* + 0j)] X= Max X 
J = 1 i 
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DUALITY RELATIONS 


A study of the geometric structure of the primal and dual models shows the 
constraint regions for both to be polyhedral. The failure modes of the frame are 
given by the extreme points of the kinematic constraint region and the facets of 
the static constraint region. Duality transformations of LP actually map extreme 
points of one model to the hyperplanes of the other and vice versa. More 
generally, 1-dimensional subspaces of the model in R n are linked to (n-1) 
dimensional subspaces of the other model. 

For primal, the random material properties specified by vector Mp occur only in 
the objective function and these determine the failure mode of the frame since 
the solution must surely belong to at least one of the extreme points. 

Therefore, unlike other procedures that require repeated solutions of structural 
model, one just needs to explore the polyhedral region for identifying the 
solution for a different (random) Mp vectors. Similarly, the vector Mp occurs 
only in the right hand side of static LP model and failure modes (facets) can be 
generated from the dual variables . 


■ DUALITY 

Static/Kinematic 

Mathematical programming 

■ PHYSICAL INTERPRETATION 

Limit states/Failure modes 

Hyperplanes/Extreme points 
Relationships 
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RIGID PLASTIC FRAMES WITH MULT I PARAMETRIC LOADING 


Proportional loading indicates a system of concentrated loads, each of which is 
proportional to a parameter, X . However, the actual loading on the structures 
may not satisfy the restriction of proportional loading. It is necessary, in 
such cases, to consider the independent variation of load factor parameters for 
the various loads acting on the frame. A static raultiobjective linear 

programming (MOLP) model has recently been formulated. Q (q — 1, 2, Q) 

denote independent load parameters and C^j and D^q are the constant coefficients. 
Unlike scalar optimization problems, the vector optimal solutions are not 
completely ordered and there is no unique ’optimal’ solution. The notion of an 
optimal solution is replaced by the concept of weak noninferior solution. 

The geometrical structure of the MOLP static model shows that it has two 
different associated polyhedra instead of just one, as in LP models. These 
polyhedral regions are defined by the feasible regions of the MOLP model in the 
objective (load) space and decision (basic variable) space, respectively. The 
two polyhedral feasible regions have frontiers made up of only polyhedral facets. 
It has been shown [6] that maximal facets of these polyhedra correspond to the 
failure modes of the structure and union of all maximal facets gives the global 
limit surface for the frame. 


STATIC MOLP MODEL 


Max 


X = Max { X , 



} 


T 


Q \ 

2 D kq q = 0 

q = 1 

k = 1, • ■ M 
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RESPONSE REGIONS 


These response patterns for frames define the limits of variation of random 
variables and once such patterns are generated for a given structure, different 
solutions corresponding to any random vector must belong to the defined regions. 
Therefore, by considering powerful algorithmic methods developed in mathematical 
and computer science literature for extremum problems, alternative structural 
responses can be predicted without reanalysis of the structure. Often, it is 
possible to further simplify the procedure in some cases, based on the 
decomposition of parametric space. These procedures decompose the parametric 
space into mutually exclusive (non-overlapping) and collectively exhaustive 
subdomains corresponding to various failure modes [7], This enables one to 
replace the consideration of an infinite number of parameter combinations with a 
finite number of parametric regions, which are also polyhedral. Multiparametric 
procedures lead to partitioning of both the load and basic variable space. All 
these procedures do not in any way depend on the probabilistic information. The 
probability considerations can be subsequently introduced to evaluate structural 
system reliability. This facilitates investigation of different loading 
conditions and probabilistic assumptions since reliability evaluations can be 
obtained without any further structural analyses. 
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SYSTEM RELIABILITY EVALUATION 


System reliability evaluation of frames for ultimate collapse is simplified by 
the use of structural responses generated by extremum procedures. For example, 
in the static case, a method has been proposed [8] which replaces the safe 
polyhedral response region of the frame with an analytically tractable region of 
equivalent volume, where the term volume is to be interpreted in a broad sense, 
since the volume may be of different dimensions and order. Reliability can be 
computed from the properties of the substituted region, which can be a 
parallelotope , hypersphere, hyperellipsoid, or any other suitable form. Use of 
hyperspherical equivalent region leads to the expression for structural system 
reliability in terms of the chi-square distribution. 

System reliability evaluation of frames for ultimate collapse by the kinematic 
approach requires the enumeration of the failure modes, calculation of the 
probability of failure for each mode and then computation of the overall 
reliability by suitable aggregation. A simulation approach that exploits the 
special structure of the kinematic model has been proposed [9] . Load and 
resistance proportionalities are determined by each simulation, and the 
associated failure mode is identified as an extreme point of the LP model. The 
procedure gains its efficiency from the fact that every simulation derives an 
associated failure condition and its probability which are then combined into a 
system reliability. 


■ SYSTEM RELIABILITY - STATIC APPROACH 

Replacement of response regions by an 
equivalent region 

Concept of equivalence 

Hypersphere, parallelotope, hyperellipsoid 

■ SYSTEM RELIABILITY - KINEMATIC APPROACH 

Generation of failure modes 

Simulation procedure 
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CONCLUSIONS 


Extremum methods of structural analysis offer significant promise for advances in 
the analysis of random structural systems and their reliability assessment. The 
complexity of the physical problem and the randomness of the variables makes the 
solution otherwise intractable. Fortunately, the mathematical nature of the 
problem lends itself to mathematical programming formulations and use of powerful 
algorithmic procedures. This has been illustrated by consideration of rigid 
plastic frames subject to collapse by flexural action. Linear and multiobjective 
linear programming models were discussed for structural systems analysis under 
proportional and multiparametric loading, respectively. Duality relations 
between the static and kinematic approaches for each of these models and their 
response patterns lead naturally to alternative methods for system reliability 
evaluation. J 


Author’s ongoing research aims to demonstrate the use of extremum methods for the 
reliability analysis of different structural systems for varying material 
behavior, structural dynamics problems and stability analysis. It is planned to 
explore structural behavior patterns with the objective of gaining insights into 
random structural behavior, dual relationships of patterns from static and 
kinematic considerations, causes of redundancy and the feasibility of using 
insights for the development of simplified and efficient computational methods 
for structural reanalysis and system reliability evaluation. 


■ PROMISE OF EXTREMUM METHODS 

■ EXTENSIONS 

Other structures 

Material behavior models 
Structural dynamics problems 
Stability analysis 

■ REDUNDANCY 

■ STRUCTURAL REANALYSIS 

■ SIMPLIFIED METHODS FOR SYSTEM RELIABILITY 
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Introduction 

A quantitative methodology is being developed at JPL for assessment of risk of failure of 
solid rocket motors. This probabilistic methodology employs best available engineering 
models and available information in a stochastic framework. The framework accounts for 
incomplete knowledge of governing parameters, intrinsic variability, and failure model 
specification error. Earlier case studies have been conducted on several failure modes of 
the Space Shuttle Main Engine (refs. 1,2,3). This paper describes work in progress on 
application of this probabilistic approach to large solid rocket boosters such as the 
Advanced Solid Rocket Motor for the Space Shuttle. Failure due to debonding has been 
selected as the first case study for large solid rocket motors (SRMs) since it accounts for a 
significant number of historical SRM failures. Impact of incomplete knowledge of 
governing parameters and failure model specification errors is expected to be important. 

Debond Failure in SRMs . 

SRM failure modes generally fall into the categories of debonding, nozzle failure, 
propellant cracking, combustion instability, field joint failure, O-ring failure, and case 
burst. As an initial case study, this work is focussing on failure due to debonding. 
Motivation for looking at debond failure is clear, as stated in reference 4: 

It is probably a conservative estimate that well over half of all mishaps (and this 
includes the latest space shuttle disaster) are due to the flame front prematurely 
reaching the chamber walls, or getting into places where it should not Usually 
the cause is a failure of the propellant-liner bond, or the propellant-to-insulation 
bond, and sometimes insulation to chamber wall bond. 

'Hie problem of solid propellant debonding has received considerable attention in the 
literature. For example, in reference 5, a finite element computer code was developed for 
evaluation of the state of stress in solid propellant case liner bond regions. Also, a closed 
form fracture mechanics solution which accounts for the dissimilar material properties on 
either side of the bondline was proposed for predicting debond growth. Reference 6 
discussed both stress-strain and fracture mechanics techniques for predicting bondline 
failures. 1 ° 

The underlying chemistry and environment of the bond region are quite complex. These 
issues are addressed in references 7 and 8. However it may not be necessary to incorporate 
all of the complexity considered in these references in order to satisfactorily assess 
reliability. Inherent variability of various parameters in the bonded region may be 
accounted for through the statistical approach described briefly below, and in detail in 
reference 1. 

The Debond Failure Mode . The sequence of events leading to failure by debonding is 
shown in Figure 1. Failure due to debonding begins as a defect, or crack, in the bond 
region, referred to as the "Initial State" in Figure 1. This defect may occur as a result of 
normal manufacturing processes, or perhaps as a result of foreign particle inclusion. In 
general there will be a distribution of size and locations of a number of defects. These 
defects may grow prior to launch as a result of induced bondline stresses from shipping 
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and handling loads, thermal cycling, and moisture absorption. Additional defect growth 
may result from vibrations during launch, axial acceleration, case 5 

aerodynamic heating, or vibrational loads. The result can be adefectof a certain size and 
location such that the flame front can enter the defect region. The debonded re S 1 ® 1 ? , 

contributes additional surface area for burning. This can lead to uneven bum and increased 
pressure. If the pressure becomes higher than the design pressure, it c ^ c ^, s ®^^ n C 
deformation and further defect growth. Case bum-through becomes a possibility, an 
detonation can result if the pressure rise is rapid enough (ret. v). 

The Probabilistic Failure Assessm ent (PFA) Methodology- 

The PFA Methodology developed at JPL is a quantitative technique for estimating reliability 
warranted by the available information. (See reference > 1 for a detailed des^ptionof PFA.) 
For cases of unacceptable risk, PFA identifies areas where design improvement and/or 
additional data are required. 

The core of the PFA approach consists of analytical engineering models which characterize 
failure phenomena in terms of governing parameters. Such J^ els f 

a failure parameter such as burst pressure, flaw size, or flaw gro**^^ 

"drivers." These drivers, i.e., the governing parameters, detenmne the value of . 

parameter. The drivers usually include geometry and dimensions, loads and environme 
conditions, and relevant material properties for the operating environment. 

In this probabilistic approach, the drivers are characterized byprobabilitydistributions 
These probability distributions express uncertainty regarding driver values within die 
ranges of possible values. The accuracy of the failure model is treated as another driver 
which is probabilistically characterized. The probability distributions for ^ <frivere 
derived from available information regarding uncertainty of their values. Thc^ fivers are 
characterized using the information that exists at the time of the analysis. There 
specific information requirement for any driver. 

The driver distributions reflect incomplete knowledge and limited 
driver values as well as intrinsic variability. The criteria of not overstati g 
information in the driver probability distributions must be observed in order to 
appropriately represent the risk that results from mcomplete knowledge and limited 

information. 

Performance, weight, and cost requirements that propulsion systems must meet may not 
permit consistently, verifiably conservative values for amdyas P^2nbv°n^^f 
cases. Deterministic analyses for such applications must be ^^ted by ot 
directly relevant past experience with applications that are similar in terms of the knowledg 
of input parameters, the validity of engineering models under the conditions of an 
application, and variability of manufacturing processes. 

When deterministic analyses are used in applications that are removed from greedy 
relevant experience base, as is often the case for launch vehicle propulsion systems, the 
uncertainty^ ur risk associated with their results increases. Since a deterministic analysis 
pmvMiT,uan.i,atfve risk estimate, an assessment of risk M as a . «suh of havrng 
chosen any specific set of values for the governing parameters of the analysis must be lett 
to the vicissitudes of judgment formed in the absence of directly relevant experience. 

Deterministic analyses of failure modes under such conditions of . U ^^i5™^ thUS 

becomes arbitrary. Launch vehicle propulsion systems are invariably 

number of failure modes for which the governing parameters may not be well known, e.g.. 
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the knowledge of loads and/or local environments may be significantly uncertain and the 
validity of analytical models used to characterize failure phenomena may be questionable. 
Under such conditions of limited information and uncertain analyses, the implicit 
consideration of risk by means of qualitative judgments based on deterministic analysis of 
failure modes is inadequate. 

In contrast, the PFA Methodology quantitatively accounts for driver specification error 
through appropriate formulation of the driver probability distribution. Application of a 
Monte Carlo technique using the driver distributions, coupled with the engineering model 
produces a set of simulated failures. These simulated failures are then fit to a parametric 
failure distribution which is treated as a Bayesian prior distribution. This prior distribution 
is then modified using Bayesian updating to incoiporate test and flight experience. The 
result is a posterior probability distribution for the failure mode. Overall mission risk can 
be estimated by aggregation of critical failure mode results. 

The resultant risk may be judged to be acceptable or unacceptable. If risk is unacceptable, 
the framework of the PFA analysis facilitates the procedure for choosing actions which will 
reduce the risk. The effect on risk, for example, of acquiring additional data, improving 
the engineering model, or making design changes, can be determined directly and 
quantitatively. 

Applica tion of PFA Methodology to Debondinp. 


At this writing, a flowchart for the engineering model for debonding has been developed. 

1 he model, shown in Figure 2, incorporates the processes described in the debond failure 
mode description above. Standard nomenclature is used: Ki and Kn are the mode I and 
mode II stress intensity factors, respectively. It is expected that some parts of the flowchart 
will require more detail while other parts represent unnecessary detail, and will be revised, 
ror example, finite element and finite difference calculations are incorporated in the 
flowchart loop. TTiese would be extremely demanding of computational time if they were 
required to be within the Monte Carlo analysis. Considerations of this type have been 
encountered before (refs. 2 & 3), and some techniques have been found which help 
minimize epu time. Modifications of the Monte Carlo approach and alternative methods 
will also be considered. 

It is important to reiterate that the risk assessment will be made using available information- 
-no additional program to develop data is required (although advantage will be taken of 
such opportunities, in particular appropriate information from the Solid Propellant Integrity 
Program, ref. 10, will be utilized). For example, variability or scarcity of data in material 
properties can be accounted for statistically. If the resultant risk is unacceptable, the 
structure of the PFA methodology can suggest options which will have the greatest impact 
on the nsk estimate. Possible options include design or processing changes, improved 
inspection capability, acquisition of additional materials characterization data, and reduction 
m uncertainty of engineering models. 


Interaction with experts in the SRM industry have been and will continue to play an 
important role in technical development of this program. 
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DEBOND FAILURE PROCESSES 


• related to inspection capability 

• spatial and size distribution of defects 

• size/geometry of grain 

• material properties 

variability: spatial, sample -to- sample 

• "stress-relieving” insulation 


* thermal environment 

• frequency, amplitude, direction of loads 


\ 

— — | • frequency, amplitude, direction of loads 

LAUNCH LOAD ENVIRONMENT * single/dual grain 

prior to combustion in defect 


• chamber pressure 

• (erosive) burn rate 

• grain geometry 

• crack geometry 

• dual grain 


• case geometry 

• case material properties 

• pressure and temperature in crack 


Figure 1. Diagram of debonding processes. Comments to the right of the boxes represent 
a partial list of relevant factors. 
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Figure 2. Flowchart of PFA debond model. 
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PROBLEM DEFINITION 


The problem is to find a global minimum for the Problem P. Necessary and 
sufficient conditions are available for local optimality. However, global solution can be 
assured only under the assumption of convexity of the problem. If the constraint set S is 
compact and the cost function is continuous on it, existence of a global minimum is 
guaranteed. However, in view of the fact that no global optimality conditions are available, 
a global solution can be found only by an exhaustive search to satisfy Inequality (5). The 
exhaustive search can be organized in such a way that the entire design space need not be 
searched for the solution. This way the computational burden is reduced somewhat. 


Problem P: Find a design variable vector x to minimize a 
cost function 

f(x) for x £ S C R n (1) 

where S is the constraint set defined as 
S - (x I gj(x) = 0, i = 1 to p; gj(x) < 0, i = (p+1) to m) (2) 

Local Minimum x* 


f(x*) < f(x) for all x 6 N(x*,8) n S 

(3) 

N(x*,5) = {x 1 llx-xll < 8} 

(4) 

Global Minimum x* 


f(x*) < f(x) for all x e S 

(5) 
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GLOBAL OPTIMIZATION ALGORITHMS 

Most global optimization methods developed in the literature are for the 
unconstrained problems. It is generally assumed that the constraints can be handled by 
adding a penalty term to the cost function. Therefore, unconstrained algorithms can be 
useful. Some of the methods date as far back as 1960s. In the following, we outline some 
of the algorithms that have been presented in the literature. 

The Tunneling Method 

The tunneling method was initially developed for unconstrained problems and then 
extended to the constrained problems [1]. The basic idea of the method is to first find a 
local minimum x* for the function f(x). Any reliable and efficient method can be used in 
this step. Once this has been done another starting point is found that is different from x 
but has cost function value as f(x*). This can be expressed as a problem of finding root of 
the nonlinear equation 

f(x) = f(x*) (6) 

that is different from x*. Again, any reliable and efficient algorithm for finding roots of 
nonlinear equations, such as the stabilized Newton's method can be used. Once the root of 
Eq. (6) has been obtained, the method to determine local minimum of f(x) is used to 
determine the new solution point. The process is repeated until there is no other root of Eq. 
(6) except x = x*. The nonlinear function defined in Eq. (6) or its modification is called 
the tunneling function. The phase where root of Eq. (6) is sought is called the tunneling 
phase. 



GLOBAL OPTIMIZATION ALGORITHMS (Cont'd) 
Stochastic Methods 


These methods are based on statistical concepts [2,3], 

Random Search 

Global-Local Phase 
Multistart 

Region of Attraction 
Clustering Method 


The Annealing Algorithm 

This algorithm, also based on probabilistic concepts, can be used to find elobal 
optimum solution [4], 6 

Can be used for discrete variables 

Analogy between Combinatorial Optimization and Annealing Process 
Concept of Statistical Mechanical System 
H(x): Hamiltonian (total energy) 

Boltzmann-Gibbs Distribution: 

p(x) = ^ exp{-H(x)/T} 


where T is a temperature and Z is a normalization constant (statistical sum). 


Let x* 


be the equilibrium configuration of the system. 


H(x*) 


mm 

xeS H(x) 


i.e.. 


Then the probability of the equilibrium state is maximal, i.e.. 



max 

X£S P (X > 


( 8 ) 

( 9 ) 


The Genetic Algorithm 

This method is also in the category of stochastic search method, such as the 
simulated annealing [5,61, in that both methods have their basis in natural processes. 

Suitable for Discrete Variable Optimization 


Three Operators: 
L Reproduction 

2. Crossover 

3. Mutation 
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ZOOMING ALGORITHM FOR GLOBAL MINIMUM 

SOLUTION 

This new global minimization algorithm combines a local minimization algorithm 
with successive refinements of the feasible region to eliminate regions of local minimum 
points to "zoom-in" on the global solution. The basic idea is to initiate the search for a local 
minimum from any point - feasible or infeasible point. Once a local minimum point has 
been found, the problem is redefined such that the current solution is eliminated from any 
further search. The search process is reinitiated and a new minimum point is found. The 
process is continued until no other minimum point can be found. 


Once a local minimum point has been obtained, the problem is redefined by adding 


an additional constraint as follows: 


minimize f(x) 
subject to 

gj(x) = 0, i = 1 to p 


gj(x) < 0, i = (p+ 1 ) to m 


f(x) <yf(x*) 

(10) 

where f(x*) is the cost function value at the current minimum 
point and y is any number between 0 and 1 if f(x*) > 0, and y > 
1 if f(x*) < 0. Constraint of Eq. (10) can be written differently as 
follows: 

f(x) < C 

f(x) < f(x*) - r lf(x*)l 

(12) 

where c < f(x*) and 0 < r < 1 . 



EXAMPLE ILLUSTRATING THE CONCEPT 

Minimize f(x) = -(x r 1.5) 2 - (x 2 - 1.5) 2 
subject to 

x l+ x 2 - 2 < 0 

-X| < 0, -x 2 < 0 

There are three local minimum points: 

1. (0,2), f = -2.5 

2. (2,0), f = -2.5 

3. (0,0), f = -4.5 

The figure illustrates the basic concept of zooming algorithm. 



Figure: Graphical Solution for the Example 
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NUMERICAL EXAMPLE 


Minimize f(x) = 9x^ + 18xjX2 + 13x2 
subject to 

2 2 

X 1 + x 2 + ^ x l = ^ 

This problem has two local minimum points: 

1. (2.5945,-2.0198), f= 19.291 

2. (-3.7322,3.0879), f = 41.877 



Figure: Graphical Solution for the Example 
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NUMERICAL EXAMPLE 


3 2 

Minimize f = 2xj + 3x 2 - Xj - 2x 2 
subject to 

X 1 + ^ x 2 - 6 
5xj + 2 x 2 < 10 

xj, x 2 > 0 

This problem has four local minimum points: 

1. (0,0), f=0.0 

2. (2,0), f = -4.0 

3. (0,2), f = -2.0 

4. (1.38462,1.53846), f = -0.003654 



Figure: Graphical Solution for the Example 





SIMULTANEOUS CONTROL AND DESIGN OF STRUCTURES 
Problem Formulation [7,8] 

State Equation: x = Ax + Bf 


0 

I 1 

R — 

0 1 

> 0 2 1 


9 B “ 

(j) T D 


Performance Index: 

PI = Jo [( X >Q X ) + (f,Rf)]dt 

State Feedback Control Law: 

f = -Gx, G = R 4 B T P 

a t p t - pbr'b t p + PA + Q = 0 

Close-Loop System: X = Ax 

A = A-BG 


Complex Eigenvalues and Damping of Close-Loop System 


EXAMPLE: ACOSS-IV Model 


Minimize weight, W = X PjAjLj 
subject to < 0, j = 1,2,.. 

co| - 5j < 0, j = 1 , 2 ,.. 

A a - Aj < 0, 

“1 = 1-341, w 2 >1.5, = 0.1093, i = 1 to 4 

For Global Solution: W < W* 


1 



Figure: ACOSS-IV Model 
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RESULTS FOR 12-BAR ACOSS-IV MODEL 


Problem No. — > 

1 

2 

3 

4 

5 

Cost Constraint 

rw*i 

100.0 

28.00 

24.00 

20.75 

19.00 

V T J. 

Optimum Weight 

31.25 

28.00 

23.29 

20.75 

No sol. 

No. of Iterations 

35 

26 

36 

28 

35 


Starting Point for all problems: 

Aj = 1000 for i = 1,2, 5,6; Aj = 100 for others 

Aj^ = 10, i = 1 to 12 
Convergence criteria: 

Constraint Feasibility <0.1% 


IlSearch Directionll < 0.01 




CONCLUSIONS 

1 . Zooming algorithm for global optimizations appears to be a 
good alternative to stochastic methods. More testing is 
needed. 

2. A general, robust, and efficient local minimizer is required. 
IDESIGN [9] was used in all numerical calculations which is 
based on a sequential quadratic programming algorithm. 

3. Since feasible set keeps on shrinking, a good algorithm to 
find an initial feasible point is required. Such algorithms 
need to be developed and evaluated. 
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FILMED 


Abstract 


One-point and two-point exponential functions have been developed and proved to be verv effective 
approximations of structural response. The exponential has been compared to the linear, reciprocal and 
quadratic fit methods. Four test problems in structural analysis have been selected The use of such 
approximations is attractive in structural optimization to reduce the numbers of exact analyses which 
involve computationally expensive finite element analysis. 


1. INTRODUCTION 

The use of detailed computationally expensive, finite element models has motivated researchers to 
ve op approximations of structural response. These approximations are useful for re-design particularly 
with use of optimization techniques, where the number of finite element analyses can be significantly 
reduced. The problem considered here is to construct local approximations using function values and 
derivatives of the structural response at one or two design points. The term local approximation used here 
!n &t a PP roxima J°n K -valid only in the vicinity of the current design point and is different from 

global approximation methods based on simplified design models or reduced basis techniques which seek to 
approximate the response in the entire design space. 

vecmr ^ f "? CUTOM I*"™. where *-<*,. *2 x„)T is a design variable 

vector Let g (x) be a structural response such as element stress or fundamental frequency, which enters as 
a constraint function m an optimal design foimulation. The problem is to construct a local approximaS 
ga(x), based on the function value and derivatives evaluated at x° and possibly another design point Then 

n S T e "?r dUatl0nS ° f C u e s u truc |, uraI response in the neighborhood of x° can be estimated using g a rather 
than the exact response g which will involve finite element computations. A variation of this problem is as 

follows: Let p be a direction vector in design space which has been determined to be desirable in terms of 
reducing the cost function subject to constraints. Usually, p is determined by solving a linear program or 
quadratic program in optimization algorithms. Now, let x 1 be a second design point along p such that 

||x -x I represents a move limit along p. The problem is now to develop a (local) line approximation g a (x) 
such that g a (x) = g(x) for points x along the line joining x° and x 1 , given by 


x = (1-0 x° + C x 1 . 0 < £ < 1 


( 1 ) 


V 

Here, the approximation g a is to be constructed using structural response information at x° and possibly x 

Hafik* t C f °^° n 0 f, va f 0us approximation methods has been earned out by Haftka, et. al. (ref. 1) and 
( a 2) ‘ ThC method , s include linear and quadratic Taylor series expansions involving first order and 
rSff? ° 7 ? er 5 enSlUVlty analysis (refs - 3 ‘ 5 )> approximations based on use of reciprocal design variables 
’ V c ° n Q v f xa PP r o x ™ations (ref. 8). Recently, force approximations have been used by 

appSadonl of ihe fo™ rat '°" P ° l> ' n0mials ^ fou " d Ref - > 0 ' "> P-Per. exponential 


ga 


n . 

cn x 2 

i=1 * 


( 2 ) 


™^ S !? e i red ^ com P ared . with linear > reciprocal and quadratic polynomial methods. It is noted the 
^ P™ enQal a PP roxir nauon discussed in Ref. (1) is of a different form than that in (2). The motivation for 

choosing exponential approximations of the form in (2) is discussed below. 
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2. BASIS FOR EXPONENTIAL APPROXIMATIONS 

The motivation for approximating structural response using the exponential form in (2) is discussed 
in this section, as also the basis for use of reciprocal variables and force approximations. The basis or 
most approximations comes from the equation 

° (A) = X (3 

which states that stress = element force/area. Area A is the design variable here. 

Reciprocal Variables 

The choice of reciprocal design variables is natural, since choosing x = 1/A as a variable results in 
a = c(1 / x) being linear in x: 


a(x) = P x 


(4) 


In the x-space, larger more limits can be imposed on changes in design, leading to faster convergence. 
Now, in a statically indeterminate truss, the stress function is of the form 


o(A) = 


P(A) 

A 


(5) 


The force P is no longer a constant, but dependent on design. The choice of x=l/A is still beneficial as it 
tends to linearize the stress function. In general, a first order Taylor service expansion of g(x) in the 
reciprocals of the variables yi , = 1 /xj, i=l, ... , n, written in terms of the original variables, xj, is given by 


g a (x) = g(x°) + X (xj -x io ) 



dg / dxj 


( 6 ) 


Force Approximations 

The idea here is to approximate P(A) in (5) by Taylor series as opposed to a(A), and obtain 

dP 

P(A °) + 77 (A o) ( A ~ A o) 

°a<A) = ^ (7) 

In the case when P is a constant, the approximation yields o a = P/A which is exact. Otherwise, curvature 
information is retained in (7) and yields a superior approximation to the conventional tangent approximation 

o a = o(A 0 ) + da/d A • (A-A 0 ). 

Exponential Approximations 

The approximation introduced in this paper is now discussed. Equation (3) may be re-written as 

a (A) = P A- 1 
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Thus, the stress is seen to be exponentially related to design variable A. This is the basis for approximating 
structural response in n-dimensional space by 6 


g a (x) = c n xj* 

i=1 

where xi, x 2 , x n are non-negative design variables. Choice of constants C and a; are discussed in the 
next section. 


A second and more general basis for exponential approximations lies in the concept of 'elasticity' a 
quantity used by economists and also relevant in nonlinear stress-strain constitutive laws. Consider a 
unction g — g(x), where x > 0 is a scalar variable. The elasticity of the function is defined as 


e - d(ln S) 

8 d (In x) 


(9a) 


or. 


_ dg/g 
8 dx / x 


(9b) 


hysically, elasticity may be considered to be in the limit, the percentage change in the function due to a 
percentage change in the variable. For instance, g=x 3 has a value e g =3, and g=px _1 has e<r = - 1 The 
exponents aj in (2) may be considered to be estimates of the elasticity at the current design point. 


In this section, reciprocal and force approximation methods have been introduced using the 
fundamental equation a = P/A as a basis. Work is being done to generalize these methods to be applicable 
to frames and certain elasticity problems as well. The exponential method of approximation has both o=P/A 

?u the r C0 ^ Cep J: °[ el ^ st icity of a function. One advantage of exponential approximations of 
orm in (2) is that for C > 0, the function ga is a monomial, which opens up the possibility of geometric 
programming (ref. 1 1). r r J b 


3. CONSTRUCTION OF THE EXPONENTIAL APPROXIMATION 

The probleirus to find the coefficient C and exponents aj, i=l, ..., k, such that the approximate 

function g a (x) = C fix?' closely matches the exact function g(x) in a neighborhood of the current point x°. 
One-point and two-point approximations will now be given. 

Uoint Approximation 


. . Here ; constants C and {a;} are determined using information only at one point x°. The technique is 

r\f e ° n C * e ^ UI ? c ^ lon v ^ ue an d shapes of g a and g. This technique has been used in the context 
d ge ° metr |c Programming where general functions are reduced to posynomial form. Morris 
log^thms SCUSSeS ^ app lcatlon of thls concept t0 structural design problems. We have, upon taking 
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( 10 ) 


In g a = In 


n 

C + £ aj In Xj 
i=1 


Differentiating with respect to xj yields 


5g a /^ x j = §a 



(ID 


Note that ga(x°) = g (x°). Equating 3ga/3xj in (1 1) to the exact slope 3g/3xj at x° yields the exponents 



3g / 8xj 
g 


( 12 ) 


The coefficient C is then obtained from ga(x°) = g(x°) as 

C = g/nx-‘ |x 0 


(13) 


2-Point Approximation 

Information at two points are used to construct the exponential approximation. Let x l) be the current 
design point and x^ be a second point, which usually is a point along a desired search direction in design 
space. The quantities g(x°), Vg (x°) and g (x 1 ) are now used to determine C and {aj}. A least squares 
formulation is adopted herein. The variable C and {at} are obtained from the minimization problem 


E = — 
2 


» < 


g 0 -c nx i( ( 

< i=1 J 


gi-cnxSj 


n 

+ s 

j=i 


^-Cnx?‘*ln xj 

v a Xj w 


v 


(14) 


The minimization of the least squares objective function E is carried out using a modified Newton algorithm, 
with a Levenberg-Marquardt correction to the Hessian when descent is not obtained (ref. 12). The 
algorithm requires the gradient vector 

VE = (3E / 0C, 3E / 9a-|, ••• , 5E / 8a n ) (15) 


and the Hessian 


3 2 E/8C 2 


h e = 


symmetric 


a 2 E/ac8a r - 

a 2 E/aa?- 


a 2 E/acaa n 
a 2 E/a ai 3a n 

a 2 E / a 2 


(16) 
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below deriVatiV6S ^ COmputed from anal y ticall y deriy ed expressions. The least squares algorithm is given 

Algorithm 2-Point Exnnnpntifll 

Step 1. Choose the initial estimates of C and {aj} from (12), (13), and e 0 = 0.001 


Step 2. Solve 


(Hg + e])5 = -VE 1 


(17) 


and update C new = C + 5,, (ai) ne w = aj + S i+ i, i = 1, ... , n. 


Step 3. Evaluate E ne w If E ncw < E the set C = C ne w, (ai) = {aj} ne w, reduce e, say, e = e/10 (if e < 
e 0 , set e - £ 0 ) and go to step 2. If E new > E, then increase £ = 10.£ and go to step 2. 

foteTce! 1 Whe " rda,iVe a " d abS0 ' U,e redUCd °" S in E f0r three co " secu,ive iterations 


4 - JEST PROBLEMS AND RESULTS 

Four test problems relating to structural design have been considered. The 1 -point and 2-point 
exponential approximates developed in Section 3 are examined. Comparison of the approximation to the 
original function is done along a line joining two design points x°, x 1 , or at points x where 

x = (1-C)x°+Cx 1 (18) 

where C is scalar variable, 0 < C < 1. For comparison, the linear (tangent) approximation based on 

ga(x) = g(x 0 ) + Vg(x 0 )*(x-x°) (19 

the reciprocal-linear approximation given in (6), and the quadratic polynomial along the line given by 

g a(0 = a + b C + cC 2 (20) 

where coefficients a, b, c are obtained from g(x°), g(xl) and dg/d C (at C=0) = Vg(x«) • (xl-x°). Thus, the 
1 -point exponential, linear, and reciprocal require only g(x°), Vg(x<>), while the 2-point exponential and 
quadratic polynomial require, in addition, g(xl). The eiror between g and ga along the line is shown both 

graphically as well as quantitatively through a relative error criterion 


RELEER = Jl [(&-&)/ gif 


( 21 ) 


and a maximum error criterion 
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MAXERR = max |gj-gai| 
1 < i < N 


( 22 ) 


Above, gi = g(x (^j)) is the exact function evaluated at the i 1 * 1 discretization point along the line in (18), gm is 
the approximate function evaluated at Cl, and N, the number of discretization points, is chosen equal to 20. 

Cantilever Beam 

The axial stress function in a cantilever beam of rectangular cross section, subjected to axial and 
transverse loads, is given as 


a(x) = 


1000 

xix 2 


6000 

+ 2 
X 1 X 2 


(23) 


where xi, x 2 are the width and depth of the cross section, respectively. The choice of design points is x° 

= (1,2) T in., x 1 = (5,8) T in. 

Referring to Fig. 1 , the 2-point exponential is in excellent agreement with the original 1 unction. 1 he 
1-point exponential behaves just as well as the 2-point exponential and is not shown in the Figure. ie 
exponential approximation to o(x) in (23) is of the form 

a a (x) = 6727. 2 x f 0891 x^ 750 (24; 


The quadratic polynomial (Fig. 1), as well as the tangent and reciprocal approximations behave very poorly. 
The values of RELEER and MAXERR in (21), (22) for this problem are given in Table 1 . It is noted that 

various choices of x® and x ^ have shown the same trend. 


Table 1 . Cantilever Beam 


Approximation Method 

Relative Error 

Maximum Error 

1 -point exponential 

0.493 

13.1 

2-point exponential 

CL422 

13.0 _ 

linear (tangent) 

700.0 

0.165 E5 

Rerinroeal 

105.4 

0.022 E5 

Quadratic Polynomial 

96.9 

0.033 E5 


143 




Tension-Compression s P ri n p 


The shear stress function in 
is given by 


a spring design problem, with xj = coil diameter and 


X2 = wire diameter, 


(x)= 522^1 

n \2 


.4xi-4 x 2 J Xl 


(25) 


= ( V 0 ’ 03)T and 7 (03 ’ °* 05 ) T in - As With the beam, the 2-point and 1 -point 

anT methndf TtS r C 056 ^ lth th f original function. Table 2 provides an error summary for 

all the methods The linear, reciprocal and quadratic polynomial are poor by comparison (Fig. 2). Other 

choices of x , x show the same trend for this problem. 


Table 2. Tension-Compression Spring 


Approximation Method 

Relative Error 

Maximum Frrnr (\ 1061 

1 -point exponential 

0.091 

0 496 

2-point exponential 

0.088 

0 067 

Linear (taneentl 

2.464 

7.265 

Reciprocal 

— 1.640 

5.876 

Quadratic Polvnomial 

11.060 

3.539 


Three Bar Symmetri cal Truss 


The natural frequency of a three bar truss (ref. 3) with x 
by the function 


1 , X 2 = cross sectional areas, is described 


w(x) = 


X 1 


2a/ 2 Xl + x 2 

Two sets of design points, leading to different performances, are chosen. These sets are 


(26) 


— 0 = (3,4)T in 2 , xl = (10,5) T in 2 

x° = (5,5 ) t i n 2 f x 1 = (l,10) T in2 


SsLl: Referring to Fig. 3 and Table 3, the 2 
best approximation, with 


-point exponential based on the best fit formulation yields the 


G)(x) = 0.2877 xf- 261 xi 0342 


(28) 
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The 1 -point exponential is poorer, with 


/ \ n ocqa „0.32 0.32 

0 )(x) = 0.2635 X-) X 2 


(29) 


Set II: Along the search direction, the original function is quite flat. In fact, the 1-point exponential 
provides a relatively poor approximation because of the flat nature of the function. The quadratic po ynonua 
is best here. Even though the 2-point exponential is second-best, (Fig. 4), the best-fit nature of the 
approximation, while averaging the error, does not provide an interval within the line where the error is 
small. This may cause difficulty for designs near the optimum. Finally, the use of reciprocal variables does 
not show any advantage over the direct variables for this case. 


Table 3. Three Bar Symmetrical Truss 


A nnrnvi motion Method 

(Set 

Relative Error 

I) 

Maximum Error 

(Se 

Relative Error 

tH) 

Maximum Error 

1_Twint pvnnnp.ntial 

0.247 

0.0286 

1.386 

0.065 

9-noint exnonential 

0.059 

0.0053 

0.605 

0.028 

T inear 

0.810 

0.100 

1.210 

0.060 

Recinrocal 

0.199 

0.022 

1.902 

0.060 

Quadratic Polynomial 

0.128 

0.012 

0.074 

0.004 


Ten Bar Truss 

The ten cross sectional areas of the truss shown in Fig. 5 are the design variables. Points x° and x 1 
are chosen as the initial and optimum design obtained in Ref. (1), as 

x° = (5., 5., 5., 5., 5., 5., 5., 5., 5., 5.) T in 2 

and „ „ ' ; 

x 1 = (7.94, 0.1, 8.06, 3.94, 0.1, 0.1, 5.74, 5.57, 5.57, 0.1) T in 2 

Again both 1 -point and 2 -point exponentials provide excellent approximations. The number of 

design variables do not seem to affect their quality. Interestingly, the reciprocal approximation provides 
equally good results, but to within a certain distance from x°. Near x 1 , the reciprocal abruptly diverges 
(Fig. 6). With smaller move limits, of course, the reciprocal will be excellent for this problem. 1 he 
quadratic polynomial provides a good approximation for this problem. Error estimates are given in Table 4 
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Table 4. Ten Bar Truss 


ADDroximation Method 

Relative Error 

Maximum Frrnr 

1-Doint exDonential 

0.088 

1.528 

2-point exponential 

0.021 

0.204 

Linear (tan genti 

0.689 

7.907 

Recmrocal 

0.599 

13.980 

Quadratic Polvnomial 

0.081 

0.803 


5. CONCLUSIONS 


Exponential functions of the form C n xf* have been used to approximate structural response. Both 

1-point and 2-point approximations have been used to determine C and (aj). The 1-point involves matching 
unction and derivative values at the current design. The 2-point method is based on minimizing a least 
squares function by modified Newton’s method. The basis for exponential approximations is from two 
sources: one is from structural theory, where a = P/A can be written as a = PA-1, while the other is from 
economics, where a function g = cx^ has an elasticity equal to the exponent a. The restriction of exponential 
approximations is X[ > 0. An advantage is that the approximating function is valid for any type of structure 
or type of structural response. Further, the exponential approximations when applied to the cost and 
constraints of an optimal design problem have the potential for being used in co^uncrion^*S^c 
programming which can effectively solve the subproblem. 

functionate nmSn °*\ °n th ? f ° Ur structu . ral Problems considered have shown that the exponential 
“ S h n^ r d d e * celIent approximations, with essentially no error, even for large distances in the 
esign space. The linear, linear-reciprocal and quadratic polynomial are much inferior. 

is ess e nria n ilv n mfnHln P e,°r b1, 5 m f i nvoIvin g natiiral frequency of a 3-bar truss, it is observed that the function 
noinfhnt K^n H \ In thls ^ ase > the 2 ‘P oint exponential (based on a best-fit) does better than the 1- 
point, but the quadratic polynomial is superior. Thus, for linear or nearly linear functions the linear 

approximation ^ ^ ** P For ° ther cases > the exponential has shown to be a powerful method of 
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Figure 5. Ten Bar Truss Problem. 
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INTRODUCTION 


n engineering design practice behavior is usually predicted based on some known nominal 
design. However, when the design is fabricated it will differ from the nominal design because of 
manufacturing tolerances. In order to generate nominal designs that will still satisfy behavior con- 
straints in the presence of manufacturing tolerances, engineers resort to the use of safety factors 
over and above those introduced to account for other uncertainties (e.g. in load conditions, material 
properties, analysis modeling). The accurate selection of the values of these manufacturing toler- 
ances safety factors is dependent on the capability of the engineer to determine the sensitivity of the 
critical constraints to changes in the design variables. This process usually leads to overly 
conservative designs. 

The task of choosing safety factors is much more difficult in structural synthesis because- 1) it 
is not known which constraints will be active at the final design, 2) as the design changes during the 
synthesis process the sensitivities of the constraints with respect to the design variables also change, 

° f thC SafCty faCt0rS themselves may change the set of critical constraints, 
ese difficulties can be overcome with the approximation concepts approach to structural synthesis 
by buffering the approximate constraints with quantities that are related to the design variable toler- 
ances and the accurate sensitivities of the constraints with respect to the design variables. Designs 
generated by this approach tend to be feasible but not overly conservative. 


Problems: 

• Design variable tolerances lead to analysis errors. 

• Choice of accurate safety factors is dependent on engineers intuition. 

Difficulties in Structural Synthesis: 

• Critical constraints in the final design are not known. 

Sensitivities of constraints with respect to design variables change during syn- 
thesis process. 

• Safety factors may change the set of critical constraints. 

Solution: 

Use the approximation concepts approach to structural synthesis with the con- 
straints buffered by values that are related to the constraint sensitivities and 
design variable tolerances. 


Figure 1 



MATHEMATICAL PROBLEM STATEMENT 

The structural synthesis problem is stated as: Minimize the weight of a structure (WO that is a 
function of the design variables (Y) subject to m constraints (displacements, stresses, and frequen- 
cies) (gj{ Y)). The n design variables are member cross sectional dimensions and nonstructural 
masses. The design variables are constrained to be in some specified range (T, L < T, < Y, ). 

The design variables have tolerances ±AT,. These tolerances may be a percentage (k t ) of the 
current design variable values. 

Minimize W(Y) 

gj( Y) <0 j = 1,2,3, ...m 

subject to 

Y\ < Y t < Y\ i = 1,2,3,.../? 

with design variable tolerances ± 

which may be a percentage of the current design variable value: 

± AY, = ± kJi 

Figure 2 
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APPROXIMATION CONCEPTS APPROACH TO 
STRUCTURAL SYNTHESIS 


nmh . n 6 a PP rox ™>°on concepls approach to structural synthesis an approximate optimization 
problem is constructed and solved at each design iteration. The use of approximations that better 

Ta^el 4 aV ‘° r aC ' Ual Pr0Wem WUI ' in genCral lead to faster desi «" convergence. Linear 
s Pproximations were first used to form approximate problems fRef 1) It was 

‘ hat displa “ men,s “ d s,resscs functions of the reciprocals of the design variables in 

Te del e l," a m rc i 11118 led 10 use °* approximations with rcspect tole inverse of 
the design vanables (Ref. 2). The use of a mixture of linear and reciprocal (hybrid) approximations 

TUndTdf SI8 ” ° f 1 Paltial derivatives ’ was found 10 •» a more conservative approximation (Ref. 
3) and led to a convex design space (Ref. 4). 

earitie!tfThe° mP hr ^ “T™ n0nlinear approximations, which capture certain explicit nonlin- 
eanties of the problem, can be constructed if approximations arc formed with respect to intermediate 

(Swir? S “ Ch SeCti °" pr0perties < Ref - 5 >’ and if intermediate response 

q tides (Ref 2) such as member forces for stress constraints (Ref. 6) and modal energies for fre- 
quency constraints (Ref. 7) are approximated. g 


Linear Approximation: 

J?t(Y) = #(Y.)+ t ^(Y,-Y ol ) 

(I) 

Reciprocal Approximation: 

*«(V) = «<Y.)+t 

(2) 

Hybrid Approximation: 

Mixture of Linear and Reciprocal Approximations based on 
the sign of ~ (assuming Y, > 0) 


Intermediate Variables: 

AV(V) = g(Y a ) + I^(X y (Y) - Y/YJ1 

(3) 

(Unear Approximation) 
For Beam Bending: 

Hiy(h, b) = g(h 0 , b 0 ) + £ ) - , o] ) 

(4) 

Intermediate Response Quantities 

««(V)=^-1.0=^-1.0 

(5) 

(Forces in Beam Bending) 

where M(Y) = M a ( Y) + £ Z22(y. - yj 

(6) 


Figure 3 
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FIRST ORDER CONSTRAINT BUFFERING 


One approach used to buffer the constraints, introduced in Ref. 8, is to add a padding term to 
the constraint function that is equal to the sum of the absolute value of the tolerance on each of the 
design variables multiplied by the sensitivity of the constraint with respect to the design variable 
(Eq. 7). This approach has the advantage that it gives good results when the constraints are nearly 
linear in the design variables. The drawback to this approach is that when the constraint function is 
nonlinear, due to the use of intermediate design variable or intermediate response quantity concepts, 
the padding term is still a linear function of the tolerances on the design variables. This can lead to 
designs that are not conservative enough. 

Another drawback to this approach is that the first order derivatives of the constraint functions 
may contain second order quantities if intermediate design variables and response quantities are 
used. These second order terms cannot be neglected since they can be larger than the first order 
terms. Calculation of the second order terms can be quite difficult, since the analytical expression 
can be very complex. The finite difference technique can be used to calculate the second order 
terms; however the error associated with this technique may become large, especially if the first 
order derivatives where generated by finite difference. The second order terms could be approxi- 
mated by using an approximate Hessian matrix (see Ref. 9), but there are also errors associated with 

this technique. 


<T(Y) 

= S(Y)+X 
/ = 1 

A Y* (Y) 

dr, 

m) 

= S(Y„)+X 
1 - 1 

3.?(Y) 
dY , 1 ‘ 

m ) 

= *™+%- r 'A¥ 


i = l 


A Y. 


cfcOO 

' dYi 


(i/r, - l/rj + X 

l = 1 


( v Y 


AO/r.) 




%00 

ar.- 


( 7 ) 

( 8 ) 
( 9 ) 


Figure 4 
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BUFFERING OF NONLINEAR CONSTRAINTS 


A more accurate buffered constraint, which captures the explicit nonlinearity in high quality 
approximations, can be constructed by using the values of the design variables at their upper or 
lower tolerance values, depending on the sign of the derivative of the constraint with respect to each 
design variable, in the constraint function. For example, in a structure with displacement constraints 
t e values of the member cross sectional dimensions at their lower tolerances would be used in the 
constraint function. The lower tolerance value is used, as opposed to the upper tolerance, because 
the lower values lead to larger displacements (the sign of the derivative of the displacement with 
respect to the design variable is negative). Since the tolerance is included in the constraint function, 
all of the nonlinearity that is captured by the constraint function is also present in the buffered con- 
straint. Note that the accurate calculation of the first derivatives of the constraints with respect to the 
design variables is the same and as simple as the method that is used for unbuffered constraints. The 

only difference is that the value of the design variable is replaced by its upper or lower tolerance 
value. 

If constraints are formed using intermediate design variables, then the values of the design 
variables at their upper or lower tolerance values are used to calculate the buffered value of the inter- 
mediate design variables. Note that in some cases, such as frequency constraints, some of the design 
variables associated with an intermediate design variable, may be at their upper tolerance values 
while the others are at their lower tolerance values. 


S (Y) = g(Y fl ) 


where Y B = i 


Y, + AY, if 


L-AK, if 


jgOO 

dr, 

dff (Y) 

dr, 


> o 
< o 


*r<V) = *<Y„) + i^p (Y ?-rj 


«1(Y) = 

■ = ' dr, v 


j i_ 

y b ~ y 

* t x nt 


£/v(X(Y)) = £(X s (Y 8 )) 


( 10 ) 


( 11 ) 

( 12 ) 

(13) 


Figure 5 



EXAMPLE 


Consider a rectangular cantilevered beam of height h and width b loaded by a moment M at 
the tip. If the intermediate response quantity approach is used, then the approximate stress is calcu- 
lated using the approximate moment In statically determinate problems such as this one, this is 
trivial since the approximate moment is constant. Hence, the approximate stress is exact. 

The approximate stress is calculated using the value of b and h. The buffered approximate 
stress is calculated using the buffered values of b and h. Since the stress is greater when the values 
of b and h are smaller, the values of b and h at their lower tolerances are used in the buffered con- 
straint approximation. Note that the constraint is a nonlinear function of the design variable toler- 
ance (Eq. 18). The buffered value of the stress is also exact. Therefore, when the design is 
fabricated and the manufacturing tolerances are at their lower values, the stress constraint will not be 
violated. Equation 19 is the first order form of the buffered constraint. Although this type of buff- 
ered constraint is exact for linear approximations, there is some error when it is used with nonlinear 
approximations because the buffering is only a linear function of the design variable tolerances. 


R 

ii 

Q IQ' 
1 

9 

a = 

6 M 
bh 2 



(14) 

-B 

= 


a* = 

6M b 


(15) 

8 

> 

b B (h B f 

M b 

= M = 

M 




(16) 

b B 

< 

1 

II 

* 

h B = 

h 

-Ah 

(17) 

a fl 


6M 




(18) 

(b - Ab) (It - Ah) 2 




~P 

6M 

Ab 

f -6 M \ 

+ 

Ah 


(19) 

a 

bh 2 + 

Wh 2 ) 
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Figure 6 
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Abstract 


Paraboloid antenna surfaces suffer performance degradation due to structural deforma- 
tion. A first step in the prediction of the performance degradation is to find the best-fit 
paraboloid to the deformed surface. This paper examines the question of whether rigid 
body translations perpendicular to the axis of the paraboloid should be included in the 
search for the best-fit paraboloid. It is shown that if these translations are included the 
problem is ill-conditioned, and small structural deformation can result in large transla- 
tions of the best-fit paraboloid with respect to the original surface. The magnitude of 
these translations then requires nonlinear analysis for finding the best-fit paraboloid. On 
the other hand, if these translations are excluded, or if they are limited in magnitude, the 
eirois with respect to the restricted not-so-best-fit” paraboloid can be much greater than 
the errors with respect to the true best-fit paraboloid. 


Introduction 

Paraboloid antenna surfaces suffer performance degradation due to structural deforma- 
tion. A first step in the prediction of the performance degradation is to find the best-fit 
paraboloid to the deformed surface. The process of finding this best-fit paraboloid has 
received some attention in the past, (e.g. Refs. 1,2) but there is no clear agreement on 
a procedure that should be followed. In particular, questions that arise are whether it is 
acceptable to change the focus length in choosing the best-fit paraboloid and which rigid 
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body motions should be included in moving from the original paraboloid to the best-fit 
one. The present paper attempts to shed some light on this second question. 

The choice of rigid body modes to be considered is associated with ill-conditioning of 
the numerical process of obtaining the best-fit paraboloid. If 2 denotes the paraboloid 
axis symmetry, then the ill-conditioning is associated with translations in the X and y 
directions. Because antenna surfaces are typically shallow paraboloids, finding X and y 
translations required to move the original paraboloid to the best-fit one leads to an ill- 
conditioned set of equations. It is possible to eliminate these translations by, for example, 
setting them to be equal to the corresponding translations at the apex. However, it is not 
clear how much is lost in terms of the root mean square (rms) surface error. This paper 
shows that not including these translations can indeed result in substantial increase in 
rms errors, but that to include them one must resort to complicated and costly nonlinear 
calculations. This is demonstrated first by the simpler case of a best-fit parabola. 


Best-Fit Parabola 

The undeformed shape of the parabola is given as 

Vi = axl 


The distortions in the X\ and y\ directions are given by f and 7] (see Figure 1) so that 
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a point A moves to position A'. The best -fit parabola is given as 

V = cx2 (2a) 

c = a + b (2b) 

The two coordinate systems shown in Figure 1 have unit vectors Zj and j x associated with 
the undeformed parabola, and i and j associated with the best-fit parabola. 

The radius vector R A from the apex of the original parabola to A' is given as 

R A = (' X A + 0*1 + ( ax A + OJl (3) 

The closest point to A' on the best-fit. parabola is denoted B (Figure 1) and has the 

coordinates [x B , (a +b)x 2 (3] in the (X, y) coordinate system. The radius vector from the 
oiigin of the original parabola to j3 is given as 

Rb = Pih T- /? 2 Ji + x B i + (a + b)x 2 B j ( 4 ) 

where and /? 2 are the coordinates of the origin of the best-fit parabola. Denoting the 
angle between the a:, and X axes as we have 


z'l = z'cos/?3 - jsin/?3 (5a) 

12 = i sin p 3 + j cos (5 b) 

Using Eqs. (5) we can obtain from Eqs. (3) and (4) the error V of A' with respect to the 
best-fit parabola 


V R b R a - (x B - x A )i + [(a + b)x 2 B - y A ]j (6) 
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where (x A ,y A ), the coordinates of A! in the best-fit coordinate system are given as 

X A = (x A + 0 cos fc + ( ax\ + 1 i) sin (5 Z - fa cos fa - fo sin 
y A = — {x A + 0 sin /? 3 + ( ax\ + r/) cos /3 3 + /?i sin /? 3 — /?2 cos /? 3 

The point /? on the best-fit parabola is found by minimizing \\R B ~ R' A \\ with respect 
to Xp. In the present work this is done with a Newton-Raphson iteration using X A as an 
initial guess ( X B is the solution to a cubic equation). 

The parameters 6, / \ are found b y minimizing the root-mean-square (rms) dis- 

tortion 


(7a) 

(7b) 


i/ 2 = — f vvdx = \' Cji/ 2 (xj) (8) 

rms 2hJ.„ 

where ±h are the limits of the parabola, Cj are quadrature weight and are points where 
the deformed parabola coordinates are given. In this work the minimization was performed 
by a conjugate-gradient method using finite-difference derivatives. 

Instead of performing this nonlinear analysis it is standard practice to linearize the 
problem. First we set COS /3 3 = 1, sin#j = A* and neglect higher order terms in Eqs. 
(7) (assuming T ] , /A , /^3 are small) to get 


x A = x A +t i + ax 2 A (3$ — (3, 
y A = -x A fo + ax\ + rj — fo 


(9a) 

(9b) 
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Next we use a linear approximation to the minimum distance as follows: We set X B = 
X A in Eq. (6) and take only the component of V normal to the best-fit parabola at 
X B — %A- This normal Tl is given as 


n = 


j - 2(q + b)x A l 
+ 4(a + b) 2 x A 2 


( 10 ) 


Neglecting higher order terms the normal component of V can be written as 


v n — v-n = -v 0 -f Fa 


( 11 ) 


where 


K = (tj ~ 2ax A ()/ -Ma 2 *^ 


(12) 


and 


= [*^,-2oa:^,l,ia + 2a 2 x i A \/yJ\ + 4a 2 x 2 A 


(13) 


The rms error is now defined as 


^rraj 


Th L * 


dx 


To minimize it we differentiate Eq. (14) with respect to a to obtain 


J_ f h 

2 h J_ h Vn daj 


dx = 0 j = 


Using Eq. (11), Eq. (15) becomes 


(14) 


(15) 


Aa = f 


(16) 
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where 


A= f tt T dx = y; Cj£i(xi)£j (xi) (17a) 

J-h i= i 

/= f ujdx = ^2 CiU 0 (xi)£(xi) (17b) 

J-h i = 1 

The matrix A is almost singular so that small deformations £ , T] can result in large 
values of ft (the x-translation) and ft (the rotation). We can minimize V 2 rms with an 
additional limitation on the size of Ot of the form 

a T c < f (18) 

and this leads to a system of equations 

(A + A I)ot = / (19) 

where A is a Lagrange multiplier (chosen to satisfy Eq. (18)) and I the identity matrix. 

Best-Fit Paraboloid 

The derivations for the paraboloid parallel the derivations for the parabola given in the 
previous section.. The undeformed shape of the paraboloid is given as 

2l = ap\ (20) 

in a coordinate system shown in figure 2. The distortions in the p\,9\ and Z\ directions 
are given by £,7/ and (, respectively, so that point A in Figure 2 moves to position A'. 


165 


The best-fit paraboloid is given as 


(21a) 

(216) 


2 = Cp 2 

c = a - f 6 

The radius vector R A from the apex of the original paraboloid to A! is given as 

R'a = [{PA + 0 COS 9 A — V sin ^J*l 
+[{pA + 0 sin 9 a + V cos 0 a ] ] 1 + (ap\ + C)h 

- x> A i 1 + VaJi + z' A k 1 ( 22 ) 


The closest point to A on the best-fit paraboloid is denoted B (Figure 2) and has the 
cylindrical coordinates [p B , 0 B , (a + b)p 2 B ] in the (x, y, z) cordinate system. The radius 
vector from the origin of the original paraboloid to B is given as 


Rb — Pih + faji + Pzfa + Pb cos 0 B i + p B sin 0 B j + (a + b) p 2 B k (23) 


where now fa , fa and fa are the coordinates of the apex of the best-fit paraboloid in 
the (^i, y\ , Z\) system. The relationship between the unit vectors in the original and 


best-fit systems is given as 



where 


(24) 


T = 


1 

0 

0 


0 

COS ft 4 

sin (3. 4 


0 

—sinfa 
cos fa 


cosfa 

0 

sinfa ' 


cosp 6 

—sinP 6 

O' 

0 

1 

0 


sinP 6 

cosp 6 

0 

_ —sinP 5 

0 

cosfa _ 


0 

0 

1 


( 25 ) 
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and /? 4 , (3$ and /3q are rotations around the axes X\,y\, Z\, respectively. 

Using Eq. (25) we can obtain from Eq. (22), (23) and (24) the error V of A with respect 
to the best-fit paraboloid. 

v = Rb — $a = ( Pb cos 6 b ~ + ( Pb sin 0 B — Va)3 + ( C P% ~ *A)k (26) 


where ( Xa , VA, Z A ), 


the coordinates of 



A in the best-fit coordinate system, are given as 



(27) 


The point B on the best-fit paraboloid is found by minimizing || V || 2 with respect to 


Pb and 9 B . Doing so we obtain 


9 B — @A + <f> 


(28) 


where 


y A cos 9 a — x A sin 9 A 
tancf) _ ~r\ , — /} 

x A cos 9 A + y A sin 9 A 


(29) 


and p B is the solution of the cubic equation 


2 c 2 p% + Pb{ 1 - 2 cz A ) - (x A cos 9 b + y A sin 9 B ) = 0 


(30) 


which is closest to p A . The parameters 6, j3 u are found by minimizing the root 

mean-square (rms) distortion 


27 r 




= ^2 J J V'VpM - ^Ci^ 2 (^i,6»i) 

0 o 1=1 


(31) 
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where h is the limit of p for the paraboloid. As in the case of the parabola the minimization 
was performed by a conjugate-gradient method using finite-difference derivatives. As in 
the case of the parabola we consider also a linear analysis setting COS Pi — 1, S i Tip t = Pi 
for i = 4,5,6. The linear approximation to the minimum distance is obtained by a 
similar procedure to the two-dimensional case: We set pg = p ^ and 0 b = ^ in Eq. 
( 26 ) and take only the component of V normal to the best-fit paraboloid at p B = Pa and 
Ob — 0 The normal fi is given by 

_ _ k — 2 cp A cosO/P — 2 cp^sinO A j 

n = t f+m (32) 

Neglecting higher-order terms the normal component of V can be written as 


v n = u.n = —Vq + Pol 


(33) 


where 

"0 = (C - 2 a p A 0! p- + 4a 2 ~p A 


(34) 


and 


= [p A , ~ 2ap A cos9 A , -2ap A sin$ A , l,p A sinff A (l+2a 2 p 2 A ), - p A cos $ A ( 1 +2a V(, )| 



(35) 


and 


— [b-> Pu Pi, Pz, P\, P$\ 


(36) 
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Note that, as expected because of axial symmetry, T) and fy do not influence the error. 
The rms error is again defined as 

h 2 tt 

-i r r 

(37) 


K 


jvlp ApiB 
0 0 


and the minimization again leads to a system of linear equations. 

Aa = f 


(38) 


where 


h 2 tt 


A= f f U‘pdpde = Y^Cili(PiA)^{PiA) 

l i L = l 


0 0 

h 


/= f f Uptpdpd$ = y: CjVo(pi, 0 j)£{pi, Oj) 

{ l L = l 


(39a) 


(396) 


Results for Best-Fit Parabola 

To illustrate the problems associated with finding the vector Or which defines the best-fit 
parabola consider a distortion of the form 

£ = £c(l - cos ^) (40) 

with rj = 0. A parabola with a/h = 0.2 corresponding to focal length over diameter ratio of 
0.625 was used. The best-fit parabola was calculated for a very small disturbance t,c/h = 
0.005. The best-fit parabola corresponding to this distortion was calculated three different 

ways: 
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(i) By using conjugate gradient optimization procedure to minimize V^ ms based on the 
nonlinear expressions in Eqs. (6) and (8). The resulting surface error is denoted V NL . 


(n) By solving the linearized problem Eq. (16). The corresponding linear approximation to 

the surface distortion is denoted V^. 

(m) By solving the size-limited problem, Eq. (19) for various values of A. 

The results are summarized in Table 1. The first line shows the results obtained with 
conjugate gradient minimization of the nonlinear expression for the error. It is seen that 
the rms value of the error can be reduced by a factor of three. However, there is great 
amplification of the disturbance with the normalized translation (3\/h being equal to 
0.1571. The linear analysis based on Eq. (16) yielded similar values for the components 
of a. However, because of the large values of these components the prediction of that 
linear analysis was erroneous. The predicted rms value of 4.9 X 10~ 4 compares with a 
nonlinear value of 6.37 X 1(T 3 . Thus while the linear analysis predicted a reduction of 
the initial rms by a factor of 3 the nonlinear analysis predicted that the best-fit linear 
parabola actually Ijjmased the error by a factor of three and a half. 

The next three lines in Table 1 include size-limited solutions obtained from Eq. (19) 
with various values of A. It is seen that as A is increased the size of a decreases so that 
the linear and nonlinear predictions become close. However, this is accompanied with 
substantial increase in the best-fit rms error. The last line in the table shows a 3-variable 
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solution obtained by setting (3\ to the apex x-translation (zero here). This solution is close 
to the large- A solution from Eq. (19). 


Table 1 shows that we have a dilemma in the construction of a best-fit parabola. Linear 
analysis requires that we eliminate one of the variables or restrict the size of the solution. 
These limitations, however, substantially increase the error rms of the now ‘not-so-best-fit’ 
parabola. The alternative nonlinear solution is complex and costly. 


This type of difficulty is not encountered when the distortion does not require [ 3 \ and 
/3 3 for its correction. As an example consider a distortion of the form 


2k x 

t = ts sin — 


(41) 


The results of the nonlinear and linear solutions are shown in Table 2 for a substantial 
value of £ s //l = 0.04. It is seen that there is hardly any difference between the linear 


and nonlinear solutions. 


Results for Best-Fit Paraboloid. 

Similar results were obtained for the best-fit paraboloid for a/h = 0.2. For example, 
a distortion of the form 

( = ( s sin(7r/9/h)sin0 (42) 

was considered, and the results for ( s //l = 0.001 are summarized in Table 3. The first 
column shows the results of the nonlinear analysis coupled with the conjugate gradient 
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minimization. The rms error is reduced by about a factor of 7, however there is again 
amplification of the distortion due to the ill-conditioning of the problem with fa/h = 
0.09. The linear analysis shown in the second column produces very similar solution, 
predicting even better reduction in rms (about a factor of 10). However, when the nonlinear 
solution is analyzed using the nonlinear analysis we find that the error actually mcj^ased 
by a factor of 3. 

The next three columns in Table 3 show the size-limited solutions based on Eq. (19). It 
is seen that, as we put more and more stringent limits on the magnitude of the solution, 
the agreement between the linear and nonlinear solution improves. However, much of the 
reduction in the error is lost, so that we have a ’not-so-best-fit-paraboloid’. Similar results 
are obtained by setting (5\ and equal to X and y translation of the apex (zero for the 
example) and solving a reduced 4-variable problem. 

While this dilemma of how to calculate the best-fit-paraboid is difficult, there is a bright 
side to it. The linear analysis gives a reasonable idea of the magnitude of error reduction 
possible with the best-fit paraboloid. 


Concluding Remarks 

An investigation of alternatives for calculating the best-fit paraboloid to a deformed 
paraboloid surface was investigated. In particular we focused on the ill-conditioning as- 


172 



sociated with the translations perpendicular to the axis of the paraboloid. It was shown 
that this ill-conditioning results in disturbance amplification so that small deformation 
can result in large translations and rotations for the best-fit paraboloid. It was also found 
that eliminating the two translations or restricting their magnitude may result in large 

increases in rms errors. 

The amplification of translations and rotations for the best-fit paraboloid results m 
grossly inaccurate prediction by linear analysis of the rms error. However, the lineal- 
analysis may be less inaccurate in predicting the achievable reduction m rms error . 
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Table 1: Best-fit parabola with various fitting schemes, £ c /h = 0.005, initial error 
L'orms / h = 1.75 X 10 3 


Fitting scheme 

b/h 

0i/h 

02/h 

0 3 /h 

v L/h Vnl / h 

rms values 

4- variable 

nonlinear 

.00104 

.1571 

.00426 

.05854 


4.90X 10 -4 

linear 

0. 

.1574 

0. 

.05844 

4.90X 10 -4 

6.37X 10~ 3 

A = .00001 

0. 

.1238 

0. 

.05844 

4.90X 10~ 4 

3.90X 10“ 3 

A = .0001 

0. 

.0426 

0. 

.0146 

8.91 X 10 -4 

9.74X 10 -4 

A = .0002 

0. 

.0248 

0. 

.00785 

9.89X 10 -4 

9.95X 10~ 4 

A = .0005 

0. 

.0113 

0. 

.00267 

1.07X10 -3 

1.07X10 -3 

3 -variable 

0. 

0. 

0. 

-.00163 

1.13X 10~ 3 

1.13x 10 -3 

Table 2: Best-fit parabola with various fitting schemes, £ s /h = 0.04, 

Vorms/bl — 

8.69 X 10 3 

Fitting scheme 

b/h 

Pi /h 

02 /h 

02 /h 

VL/h v NL /h 

rms values 

nonlinear 

.0142 

0. 

-.00228 

0. 


6.02X 10 -3 

linear 

.0141 

0. 

-.00225 

0. 

5.69X 10 -3 

6.02X 10 -3 
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Table 3: Best-fit parabola with various fitting schemes, (, s /h — 0.001, initial error 


Vorms/h = 4.868 X 10 


-4 


Fitting 

scheme 


6- variable 


4-variable 


nonlinear 


linear 




O 

II 

A = 5 X 10" 

" 6 A = 2 X 10 -5 

A = 5 x 10" 

-5 

b/h 

2.33x10" 

4 0. 

0. 

0. 

0. 

0. 

A/h 

0. 

0. 

0. 

0. 

0. 

0. 

A/h 

-0.09055 

-0.09055 

-0.06281 

-0.03280 

-0.01686 

0. 

A/h 

1.41x10" 

' 3 0. 

0. 

0. 

0. 

0. 

A/h 

-0.03366 

-0.03366 

-0.02313 

-0.01174 

0.00569 

7.12X 10~ 4 

A/h 

0. 

0. 

0. 

0. 

0. 

0. 

I'/./h 


5.19x10" 

5 1.12xl0 -4 

2.14X 10 -4 

2.70X 10 -4 

3.29 xlO" 4 

tWZ./h 

6.91X10' 

_5 1.47x10" 

' 3 7.05X10~ 4 

2.79X 10 -4 

2.73X 10 -4 

3.29 X 10 -4 
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Abstract 


This paper presents approaches to the combined design of structures and controllers for achieving optimal 
maneuverability. A maneuverability index which directly reflects the minimum time required to perform a given set of 

desigmng 1,16 flexible appendages, the maneuver time of the spacecraft is minimized under 
^ ^" S “’ aintS h 0f st ™ctural properties and of the post maneuver spillover being within a specified bound. The spillover 
ach,eved by ^"8 f® of 3,1 appropriate reduced order model. The distributed parameter design problem is 
!"h as A sumed shape functions, and finite element analysis with dynamic reduction. Solution procedures have 
gated. Approximate design methods have been developed to overcome the computational difficulties Some 
new constraints on the modal frequencies of the spacecraft are introduced in the original opUmization problem to facilitate 
the solution process. It is shown that the global optimal design may be obtained by tuning the naHira*] frequencies to 
sausly specific constramts We quantify the difference between a lower bound to the solution for maneuve? time 
associated with the original problem and the estimate obtained from the modified problem, for a specified application 
requirement. Numerical examples are presented to demonstrate the capability of this approach. 

I. Introduction 

6 structure f s sucb 3S antennas or space stations will be very flexible, not only because of the high cost of 

mt f h f , H C i tUreS [ r0m 10 space ’ but 3)50 wiU he constructed or deployed in orbit and will 

r™t t0 . Wlthstand lar 8 e launching and gravity loads. However, when a space structure is very flexible its active 

fl Ste hi Ca " C . XClte and othe ™'se significantly interact with its flexible modes. Thus, the idea arises of achieving 
die best flexible mode suppression for attitude maneuver of spacecraft. The control problem of time-optimal rest to rest 
slewing of a flexible spacecraft through a large angle has been investigated [1], In that work, a specific spacecraft is 

U f ng 3 reduced ordcr model > 3,1(1 *** dme-optimal control history of this modelled system^ derived In some 
ume critical applications, it is required that the maneuver be performed as rapidly as possible. As a consequence 

r!Zh,WH °P timiZ3 | J0n 1S considered so as to further minimize the maneuver time. The whole design prwess the idea of 
combined design of controllers and structures for optimal maneuverability, is considered in this work. 

. m J 1 ? ldonal,y, f th ^ overa11 design Problem for actively controlled space structures is treated via an iterative two-part 
Kheme. Redesign of the structure including sensor and actuator placement is performed in one stage, and then the control 
aw is modified for the resulting system to complete an iteration cycle. Generally different design objectives aoolv in the 
An More r cc ent ly. the need to integrate the design of a structure and its control system has been recognized 

S5S "“K **■ ^ StrUCtUral and COntro1 ™ substantially coupled. 

Bodden and Junbns [2] presented a method for eigenvalue optimization with sequential or simultaneous design of 

3nd COntrol . Khot : 0z . Jenkayya, and Eastep [3-5] considered structural optimization, including constraints on 

SS ^T’ ^ °" 3 h"^ 33 ^ model of the controller. Ha"e, 
Lisowski, and Dahl s [6 7] treatment of the problem of simultaneous structure and control design for a maneuvering 

spacecraft resulted m a lmear-quadratic optimization problem. Bendsoe, Olhoff, and Taylor [8]presented an algorithm for 
' n ^^ deS, 8 n of jbe structure and its control system which includes a constraint to limit tire control spillover into the 
" ' ! d ? des '. Lust 3,1(1 Schmit f 9 l Presented a control-augmented structural synthesis methodology P in which the 

structural member sizes and active control system feedback gains are treated simultaneously as independent design 

, V3 " 3b *f-. ° nod33nd Haftka [10] considered the optimization of the total cost of the structure and control system subject 
* T ns ma g n ‘ ltu ^ of the response to a given disturbance involving both rigid-body and clastic modes Lim 

and Q^ kU, r 1 1 presented 311 ldea for optimizing the robustness of structures and structural controllers using homotopv 

wii Z^ r " ear programm,n 8 algorithms. Khot [12] presented algorithms for design of minimum weight strucSes 
with the goal of improving system dynamics by use of a closed-loop control system. structures 

linear °l the develo P men ‘f on simultaneous design of structures and controllers reported in the literature use simple 
linear feedback control laws and quadratic performance indices. Practical constraints such as lim£n , Stoe ZSSto 
of the control effort generally are not taken into account. Tire use of such relatively narrow f“fw^ 
may have senous implications in terms of the usefulness of the results. It is understood for example that tte uS 
jrerformance radices expressed as hnear/quadratic functionals is generally inappropriate unless loop tLisfer recovery 

2? cXri , 1 ^ r?r into "* formul3lion - Prnlhermore, the constraints usually^ literate are on 
e closed-loop eigenvalue distribution and structural frequencies. These constraints are not as direct to the application 
5255"“ “"fc °n rise time, maximal displacement, or maximal sires, 

degradation of the opumal system coming from the control and observer spillover is also generally noHncluded. 


180 



In the present work, we examine the problem of fully coupled design for a spacearaft and its associated control. 
The design of the structural system and control is to be integrated so as to optimize with respect to a single cost 
function The objective is chosen to reflect the maneuverability of this structure/control system, i.e. the time required to 
perform a given maneuver or set of maneuvers. Various forms of Mission Specification can be reflected in the definition 
of the performance index. Ours includes criteria related to sets of maneuvers with specified probability of occurrence. 
This performance index is generally more meaningful than the usual LQG index with minimum weight. The minimum 
time' objective is appropriate for application in slewing or other retargeting maneuvers. Furthermore, the problem is 
formulated in a way to accommodate in explicit form of various practical constraints, such as limits on control action and 
performance error (control spillover). Also, the formulation is consistent with a nonlinear bang-bang form of optimal 
control design. 

The spacecraft is modeled as a linear, elastic, undamped, nongyroscopic system. The necessary-and-sufficient 
condition for the time-optimal rest-to-rest control problem can be considered as a mapping from the structural dynamic 
properties to the optimal maneuver time. The maneuverability is optimized by updating design parameters. 
Characteristics of the problem and problem solving procedures have been investigated. Approximate design methods 
have been developed to overcome the computational difficulties. Numerical examples are presented to demonstrate the 
capability of those approaches. 

II. Combined Design of Structures/Controllers - 
Problem Formulation 


Consider the linearized rotational dynamics of a flexible spacecraft where control inputs are used to actively control 
the rigid body mode and flexible modes. The spacecraft is modeled as a linear, elastic, undamped, nongyroscopic system. 
There is a rigid central body, as shown in Figure 1 , to which N (N > 2) identical flexible appendages are attached with 
uniform spacing between them. Along the appendages, there might be some kinds of distributed or concentrated payload 
masses for practical usage. The spacecraft may be very large and flexible. The spacecraft is to be controlled by a single 
torque actuator located on the central body and m torquers located at identical locations on each of the N appendages. The 
amplitude of the torque applied by each torquer is limited. The objective of the control design is to time-optimal ly 
slewing the spacecraft through a specified angle 0, and achieve flexible mode suppression at the end of the maneuver. 
Assume that the appendage displacements, slopes and central body rotation rates remain small and the appendages are 
inextensible. The appendage displacements are restricted to a plane orthogonal to the central body’s axis of rotation. 


central 
hub 


slew axis 



Lump mass 
(payload) 


Design of the N 
Identical appendages 
Figure 1 


The design parameters of the appendages can be the cross section, stiffness or density of the material, layouts of 
the composite material or the location of torquer actuators along the appendage. Let the design parameter vector be E, 

€ R N , implying that the structural dynamics properties are implicit functions of 


Maneuverability Index 

The maneuverability, formulated as a maneuverability index, reflects the cost required to perform a given maneuver 
or set of maneuvers. The mission profile is specified by giving the probability density function p(9) of the required 


maneuver amplitude 0. Let tf (0) be the optimal maneuver time for maneuver 0, and t f *(0) is a function of the 
Se^ra d b%^S? rV “ tori mancavciab ^ is **> a function of §. Wedefinethe 


For example, let p(0) 5(0-0 ) then p*(£) - tf (0*). In other words, the maneuverability index represents the expected 

wiS e r^t a IoT neUVer Ume f0F 3 giVe " miSSi ° n Pr0me ' ThC StrUCtUral deSig " Pr ° blem iS thCn t0 ° ptimize 


Optimal Design Problems 


Assume that the structural design parameter § is restricted to belong to a compact set E which represents 

Assume that the design of the appendages will not change the characteristic of the torquers along the 

TW?fv? 8CS ' ^ 0th f er W ° 1 rdS ’ am P 1,tude Umits of 1,16 tor( l uers remain «t*e same for all values of the design parameters 
Therefore, we can formulate the optimal combined structure/control design problem as : 

H=mln p*(5) 

where E is the space of structural design variables 
subject to two sets of constraints : 

I. a. Material resource constraint, 

b. Geometric configuration constraints : 

such as the max. and min. thickness limits of cross section, 

c. Dynamic response constraints : 
such as the max. stress and displacement limits, 

II. The post-maneuver control spillover is within a specified bound. (2 2) 


and. 


Constraint II takes into account the performance degradation associated with the unmodelled dynamics. 

We approach the distributed parameter design of the cross section of the appendage using assumed shape functions 
For examp e let the design parameter of the cross section be the thickness d.stXtion We assumed at tte thic toess 

functions ^Thi^ ^ ? caUon ^ ong the a PP enda g e ) is represented via a linear combination of a set of assumed shape 
functions. This approach uses the same idea as design variable linkage [17]. ™ 

Tte distributed infinite degree-of-freedom system is approximated with finite elements We discretize the 
^pacecraft into a finite number of elements and then perform modal analysis. There are two kinds of mathematical 
mimhc' f °f and ana, ys is ; Let subscript E indicate a quantity derived based on the control evaluation model The 

m^rf A f m0dCS “l^ 1S ,S ^ number of degrees of freedom in the finite element analysis (let it be n in this 
paper). Assuming this model to represent the exact dynamic system, we can evaluate the performance of the controlled 

1 ^ l UteC T *'" diCa " .* -terived based on ibe conirol design mo dd. re^oltsig" is 

the model on which we obtain the optimal-time maneuver law. We assume there are r (r « n) vibrational modes 
tamed mthi s model. The natural frequency and mode shape of the modes in the control design model can be easily 
obtained from the control evaluation with the dynamic reduction method. y 


Results of the Linear Time-Optimal Control Problems 


Results presented in this section were obtained in the recent paper [1]. The optimal control character^ h™» ; c 

tetanbSuT Sy “ m ThC problem « f Ume-oplinial rest-to-rest slewing mieuvercan 

Problem M(0)^: 


min tf<0) ( 2.3) 

Subject to : 

x = A x(t) + Bu(t ) 
luj(t)l < Uj . j = 0, 1 m 

where uq is related to the control input at the central rigid body and uj, u 2 , ... u m are related to the m 
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torquer actuators along the appendages. Uj ; j = 0, 1 m are the corresponding amplitude limits. 

x (0) = (0, 0, 0) 1 , where ( ) t denotes transpose. 

x(t f ) = (0, 0, 0) 1 

A = Block diag[Aj], B = Blocked [Bj], where 


f° l ] 

Lo oj 


Bi 


0 co . 

I 

- to . 0 

I 

'0 0 . 
_p‘op; • 


0 

p'r 


= 0 , 


COi is the natural frequency, i = 1,2, ..., r 


; i = 1, 2, .... r 


(2.4) 


Let the solution for problem M(0)^be tf (0). 

* 

Theorem 1.1. Let tf* be the optimal maneuver time. For all 0, Problem M(0)^has a unique solution tf . 
Theorem 1.2. For a given 0, the optimal control law is of bang-bang type, and is symmetric around tf 12, i.e. 
u(t f */2- 1)* = - u(t f */2+ t)*, 0< t <t f * 12. 

Reference [1] treats the general multiple control case, where there are m+1 control inputs. However, for simplicity, 
herein we assume that only one control input is used to control the maneuver, that is, the scalar control case. This 
assumption means that the N torque actuators on the appendages and the actuator on the rigid central body taken together 
represent one control input. 

Theorem 1.3. Assume there are k switching times between 0 and tf /2 , and let them be tj , i = 1,2,... k. Let J* be 

the total rotational moment of the spacecraft, and (p 0. 0. P 0. 0. — » P 0> 0) be thc costate variable at mid-maneuver 
time. Then the optimal maneuver time and the switching times satisfy as necessary and sufficient conditions, the 
following system of nonlinear algebraic equations : 

(r;> 2 - 2(t k ) 2 + 2(, k J- ... +2(-l)‘<i 1 ) a -e/*/l^ (2.5) 

cos (co .r,72) - 2cos(co. r J + 2cos(co, f 2( - l)*cos(co ,.!,) + (- 1)* + ’ =0 

( 2 . 6 ) 


(2.7) 


i- 1.2, ... r 

ff/ o r*/2 U Q sin (to . t* /2) U Q tia(m r t*l2) 

o 

P 0 

P P 
K o^o 


- 1 
0 

: ► 

t sin (co,/J sin (to r r ) 

k 1 * * 

► = < 

: : i 

: 


0 

t sin (to O sin(co r /) 

ip 


.0 . 


and two inequality equations : 


t f /2>t k > t R _ 1 > >t 2 >t 1 > 0 

r 

P° 0 1 + XP'; PlSinfrit ) * 0 


i =i 


(2.8) 


where 0 <t <tf*/2 totj i = 1 . 2, ... k 

To solve for the optimal control history, we need first assume the number of switching times, say k, then try to find the 
solutions { t- , j = 1, 2, ... k J and { p j 0 , i = 0, 1, 2, ... r ) for (2.5)-(2.7). If (2.5)-(2.7) really admit solutions and 
they satisfy (2.8) as well, by uniqueness of the solution of the optimal control problem, we have the unique solution. 
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k?s Sr S’,? general : k is . always ^ 10 r - 0nly when { 0)i , i = 1, 2, ... r } satisfy some special conditions, 
k is less than r. For the case where k is equal to r, Theorem 1.3 can be simplified by omitting (2.7). 

Reduced Order Model 


.,„ r ^ C K n ° W c ° nsider h™ many flexible modes should be retained in the reduced order model. The question is 

thC dcgradaU ° n Ir \ lhc P crform ancc of the designed system on the control evaluation model This 
performance degradation is associated with unmodelled dynamics' of the uncontrolled residual modes in the control 
evaluation model (from r+l-th term to n-th term). The effects, therefore, result in post maneuver free vibration of the 

SSSEEL3? nMd 10 ^ SUr ' ' heK VibraUom ta « " l *“» * ■P*™ performance 

There are two ways to quantify the performance degradation : (i). the residual or spillover energy &(t) , and (ii). 
the pointing error of the rigid central body after completion of the maneuver 0e(t) (where t > t f */2). From the recent 
investigation on these [1], the latter is the better one because the maximum pointing error continues to decrease as we 
suppress additional modes at the final time, while the spillover energy does not necessarily decrease. Also [ 1 ] gives three 

closed form expressions for the upper bound |M'>|, based on the control evaluation model. Among them the most 
useful according to our experience is 6 ’ usl 


|e, (I l < 21 (2 + it WJJ * 2 ) £ (p‘ o > 2 ; , > ,2. 


0f ^ COm °' ra ° dd in ° ,<ter “ ° tai " a p-specificd po s , -maneuver 

Characteristics of the Optimal Design Problem 
Theorem 2. 

Suppme the number of flexible modes retained in the model is fixed. The optimal maneuver time solved from the (2 5)- 
(2.8) is a continuous function of the structural design variables, v ’ 

Corollary 1. The objective function, p(£), is a continuous function of 
Corollary 2. There exists a solution to the optimal design problem (2.2). 

rnn ci,wm haVe observed thal Ibe objective function is always a differentiable function of the structural design variables t 

S CaSC , where f k ,s ^ ual t0 r - The 0 P dm aI maneuver time can be obtained from (2.5)-(2 6) Actually ’ * 

(2.5)-(2.6) represent a system of implicit equations of the form : v y 

F (t f * , t, , (0, p, J*) = 0. (2 ]0) 

The gradient of the optimal maneuver time with respect to structural parameters can be obtained using the Implicit 

Function Theorem as follows : Let x=(tf , tj) and y=(co, p, J*) 

Theorem 3. (Implicit Function Theorem) 

Suppose (x 0 . y 0 ) is such that F(x 0 , y 0 ) = 0 and F(x 0 , y 0 ) e C k , and the Jacobian matrix [OF/dx] is nonsingular 
(regular) at (x 0 , y 0 ). Then there exist a neighborhood of y 0 , say N(y 0 ), and a mapping G : N(y 0 ) -> R n such that 
*0= G(y 0 ) and G(y 0 ) e C k , and F(G(y), y) = 0 on N(y). Moreover, we have 

[3G/9y]t | y0 = - [8Fy0x]-«| xQ [8Fy9y]1 xQ , y0 

2JS Chain Rule we can obtain the gradient of the objective function with respect to the design parameters A 
programming tcffinTif 11 ^ Kurash - Kuhn - Tu <*er necessary conditions [18]. We use mathematical 

From (2.10), t f is an implicit function of (to, p, J*). Furthermore, for the generic case where k is equal to r t f * 

actually is a function of to and J* only. We show the behavior of t f * (to, J*) for the simplest case where there is only 
one flexible mode in Figures 2 and 3. y 


184 



60.0 



30.0 


25.0 


0) 

| 20.0 

u. 

V 

> 

3 

S 15.0 
0 


0 

£ 

a 

O 


10.0 


5.0 


0.0 



Rcfiarvnca 

W1 as 1.50 


W1 = 0.50 

W1= 2.0C 

- 

W1 = LOO 





_ 




■ . i ■ ■ ■ ■ i**- 

. . i . . , . t , . . . 


0 0 0.5 1.0 1.5 2.0 

Natural freq of lat Mode (rwifrac) 

!>■ to th. nuiuw Urn. lor • 

TM Mm t j* t {QiDiunr tngk and total torqua limit uo 


0.0 20.0 40.0 60.0 80,0 100.0 120.0 140.0 160 0 
Angle*J/Tmax - flV2)**2 
Tt* is tha manw^ar Uma of rigid awwacraft with tha rotation moment 
J and tha toUl torqua Umltlteiax. 

Reference La the curve : IT ■ 2*SQRT(A ngia*J/Tma x ) 


Figure 2. Figure 3. 

We have found the following results : Assume a spacecraft has only one flexible mode. 

(i) . For a spacecarft with very small a> x (usually a very flexible spacecarft), t f * is quite large. On the other hand, for a 
spacecraft with large (as shown in the Figure 2, greater than 2.0), tf* is almost the same as that of the equivalent 

rigid spacecraft. * 

(ii) . For a spacecraft with 0 J*/Uq > 120.0 (the torquer limit is very small or the maneuver angle is very large), tf is 

almost the same as that of the equivalent rigid spacecraft. 

Of course, a typical spacecraft has more than one flexible mode, and we can not say much about it. However, Figure 2 
and 3 provide important information. If the spacecraft is very flexible or the torquer limit is very large (usually this 
implies very large maneuver speed), the result of the optimal design can provide substantial improvement. 


Problem Solving Algorithm 

The control design model is chosen according to the analysis of control spillover. In order to take advantage of 
Theorem 2, we assume that the control design model is fixed during the optimization and formulate the optimization 

procedure as 
2>i : 

Begin with a reasonable baseline design of the spacecraft. 


Step 

Step 

Step 

Step 

Step 

Step 

Step 

Step 

Step 


0 

1 

2 

3 

4 

5 

6 

7 

8 


Set up the reduced model by (2.9). (Set the value of r) 

Initialize the design variables. 

Get the cross section of the appendage for the current value of the design variables. 

Finite Element Analysis. 

Calculate the natural frequencies of the modes in the reduced model by the Dynamic Reduction 
Method. 

Solve the time-optimal control problem to obtain the optimal maneuver ume. 

Find the next values of the design variables by the Nonlinear Programming. 

If the result is convergent, Step 8. Otherwise, go to Step 2. 

If the spillover constraint (ii of (2.2) ) is satisfied, then Stop. Otherwise, Step 0. 


Although the algorithm fPi is able to solve the optimal structural design problem (2.2), unfortunately, in our 
experience, there exist a lot of numerical difficulties associated with it : 

i) . To solve the time-optimal control problem, we need to know the number of switching times. 

ii) . Actually the set of nonlinear equations (2. 5-2.7) admit many solutions, of which only one satisfies the 

inequality conditions (2.8). Thus, even though we have a good nonlinear equation solver, it would not be able to 
guarantee to find the solution we want. 


OF POOR QUALITY 





Given ail the difficulties above, it seems a formidable task to solve the optimal design problem by without any 

esp f c ‘ ally ‘ f one to find *•« global optimal design. Therefore, we introduce approximate design 

methods as described in the next section. ® 

III. Approximate Design Method 

_ Tbe fundamental idea of this solution procedure is to formulate an approximate design problem without violating 
a y C °" s | r i a , ,nt °. f P^lem. The solution of the approximate design problem is a ‘near-optimal design’ in the 

sense that there is little difference of objective function between the two solutions. We need to quantify the difference 
without solving the original problem and make it as small as possible. However, there is a trade-off between accuracy 

ull^de m/° r 8 , P S H ThuS “ ' mportant capability of the approximation algorithm is that we can adaptively 
upgrade the approximation procedure to obtain a reasonable result according to the specific application requirement 
Since design models can not exactlly represent the real system, it is unreasonable to concern oneself so much about a 
relatively small improvement of accuracy of the solution based on a design model. In this section we introduce two 
approximate design methods : the Adaptive Frequency Tuning method, and the Minorant method. The former one is 
suitable for the single maneuver case; the latter one requires more computation work but is suitable for the 
multimaneuver case. 

Frequency Tuning Approach 

There are two basic assumptions : 

;*""!* V™ } * : " atural frequencies of the modes retained in the reduced order model can be freely assigned by 

adjusung the values of the design variables. B y 

thfSlic ' hwk ,• During ih u dcsi ? n iterati °n, the mass distribution of the appendage is taken to be independent of 
e stiffness distribution, i.e the total rotauonal moment of the spacecraft, J* , does not change when the 
stiffness distribution is modified. u 


satisfy C ° nSldCnn8 (2 ‘ 5 ’ 2 ' 6) ’ ' f f ° f 3 spacecraft the natural frequencies of all modes in the reduced order model happen to 

®i * tf*= j| * 4 7C , i = 1,2 r (3_j) 

where tf is the maneuver time and j j is some integer multiplier, 
then, the solution in terms of switching times and the optimal maneuver time satisfy : 

V = n anH .* 2 /e/ *</ 


k = 0, and tf = i \' 3J v o (3 2) 

It also satisfies the inequality condition (2.8). Thus we solve the time-optimal control problem for {CO;, i = 1,2, r 

satisfying (3.1) } . Moreover, (3.2) imply that there is no switch of the control history between 0 and t f */2, and only one 

hv 'S This means 311 flexible modes in * he reduced model are dead beat at the end of maneuver 

optimi^rn problem maneuvers a r >8 id body of the same value of total rotational moment J*. We have the new 


min tf = 


/QJ *fj. 


subject to : 

the constraints I and II of (2.2), and (3. 1) (3 y 

TgtoS optimum"^' ^ ASSUmpti0nS 1 and 2 above ’ ** so,ution of ^2 solves our original problem, (2.2), and it is 

forS^ap^swte SST n “' ble m0deS VC,y Sma " X ” d Very " 8h ‘ <in lhal J * is s ™">- This idea 
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Adaptive Upgrade Algorithm 

Unfortunately, Assumption 2 above is not always satisfied in general applications. For example, in designing an 
appendage of rectangular cross section with high density material, the stiffness is highly coupled with the design of mass 
distribution. Actually, implicitly assumes that the global optimal design of the appendages is such that the time- 
optimal control is the same as the rigid-body control strategy. We restrict ourselves to solve the original problem in a 
subspace of the feasible design variable space. Therefore, the result of 2>2 in general does not apply and needs to be 
modified or upgraded. 

We first quantify the index of improvement in approximation as the difference of objective function between the 
exact optimal design (the solution of original problem) and the solution of approximate design problems. Let tf be the 

maneuver time of the exact optimal design, which is equal to the minimum of tf over the entire feasible design space. 

Also we note that (I is equal to jp (6 . ) • f/Oj) dQj . Let tf a be the approximated solution of tf and ft respectively. 

Then we introduce : Index of approximation : Tq= \ tf a - tf I or I (I - ft I (3-4) 

An approximated solution is better if the index of approximation is smaller. However, this doesn't mean the two designs 
are close to each other. For example, they may be substatially different in shape. In order to avoid difficulties in 

computing tf, we modify (3.4) : £1 = I tf - L { tf) I or \\1 - L (ft) I , 

where £ (•) is a lower bound of • , and it is very easy to compute. 

Also, we have l\\1)= Jp(9_ ) (f(9i)) d9| 

There are two ways to define such a lower bound : 

(i) . the maneuver time for a rigid spacecraft with the least fesiblc total rotational moment J*: 

It is usually unreasonable to define the lower bound in this way because (3.7) is very conservative. The appendage with 
the least total rotational moment is usually too slender, too flexible, and likely requires a long maneuver time. 

(ii) . the optimal maneuver time of the optimal design which is based on a reduced model with only one flexible mode. 
Let the superscript 1 of tf indicate that the value is based on a reduced model with only one flexible mode. Thus 

tf ) = tf 1 = minimum of tf ^ over the entire feasible design space. (3.8) 

Since we need more maneuver time for the reduced model with more flexible modes, we know tf 1 is a lower 
bound of the maneuver time for the design problem of any reduced order model. We need some computation effort to 
calculate tf h however, the calculation is not very difficult. It is more reasonable to define the lower bound of the 

maneuver time to be tf 1 . 

We propose the modified approximate problem 1^3 according to the following facts . 

Fact 1 : For a specified reduced order model with r flexible modes, we can divide the feasible design space into . 

D 0 : ( ^ ; the time-optimal control history of this design admits only one switch at mid-maneuver, without any 

switch in (0, tf* /2 ) } , urn ,^n \ 

'Di : ( ^ : the time-optimal control history of this design admits at most one switch in (U, lp/4 j, 

©2 • { £ : the time-optimal control of this design admits at most two switches in (0, tf* /2 )} , 


(3.5) 

(3.6) 

(3.7) 


and ©o — ®j — ©2 — — 


Fact 2 : tf* over D r > tf* over © r+ j 

Fact 3 : the solution of (2.2) is the tf* over © r for some r > 0. 


(3.9) 

(3.10) 


Actually, the solution of (?2 is nothing but tf over ©p. Similarly, 2*3 is the problem of solving for tf over © r , 
r > 1, adaptively upgrading with respect to the index of improvement, and with a stopping criterion based on sufficiently 
small change of improvement We can eventually obtain the exact global optimal design if the upgrade goes on. 
However, we have restricted ourselves to solveing for tf over © r , r < 2. 
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The Minorant Design Method 


!?2 and (P 3 are not suitable for the general multiple maneuver case because it is difficult to find to- , i = 1,2, ... r 

which satisfy (3.1) for many different maneuvers, {Gj}. In this section we discuss an algorithm, the minorant method, 

which is more difficult to implement, but, suitable for the multiple maneuver case. While solving the time optimal 

control problem, we find that for any design of spacecraft, t f * r+ ' > t f * r ; however, the difference becames smaller and 

^m.!i S innT reaSeS ' ^ 2 ™ ^ nu ™ e ^ stt i di< *. it is observed that the maneuverability is most influenced by the 
tot^rotauonal moment J and then from the few lowest flexible modes. An appendage with smaller total rotational 

moment or with more rigidity, in the sense that the natural frequencies of the lowest few flexible modes are large tends to 
be very maneuverable. B 


(P 4 is based on the following assumption and fact. 

Assumption : For any fesible design of the spacecraft £ e E, we have ! - tf*(£) I+i | < | tf*(£)* +i - 

t f £)» I . i > 0. where the superscript i indicate that the quantity is obtained based on a reduced model with i flexible 
modes. 

Fact 4 . I tf tf* I < I tf** 1 - tf* \ , i > 0, and I tf*- tf T I -» 0 as r and i are sufficiently large 
Furthermore, \^. p i+i | < I ^ ^ | , i > 0, and 1 4 ^ r U Oas r and i are sufficiently large. 

Step 1 : Let i = 0, and Solve (I 1 by (P <\ . 

Step 2 : Obtain the index of improvement £. If there is no relative change of improvement stop 

Otherwise, i = i+1. Go to Step 1. 

The exact optimal design can be obtained for i = r. However, we do not go beyond i > 2. The capability of Ta will be 
investigated later with numerical examples. H 


IV Numerical Examples 

■ (h exa f lpl “’ we consid er designing appendages by adjusting the cross section. We use practical examples 


In what follows, we perform the modal analysis with the finite element method, and model the flexible spacecraft 
°" e ngl d mode and twenty flexible modes. There are r flexible modes, obtained 

metiiod, retained in the reduced order model for control design. The reduced order model is specif! ed according to the post- 
maneuver asO nslT ^ examp T le K S ’ we specify thc maximum Aguiar deviation of the central ng.d bodypost 
SKS2? nfL aP f ndagCSarC l ~ bcW T (as shown in F ‘g ure 4). Our goal is to obtain the optimal flange 
i f appCndageS ’ assume the width of *© web, and thickness of the web and flange to be 
constant. The flange depth is symmetric about a central line passing through the cross section. We use two spline 
polynomials as the assumed shape functions to describe thc half flange depth : Sp " ne 


A, (Jr )= Cj + ( c 2 /L )x + (c 3 /L) 2 x 2 + 


<C 4 1L) 3 x\ 


0 <x <L/2 


L/2)+(c 5 /L ) {x 


-L/2) 2 + (c 6 /L) 3 (x 


-L/2) ,L !2<x <L 


h 2 ( x )=h i (L/2)+ h x \L 1 2){x 

where cj, i = 1,2, ... 6 are design variables. ^ 

We TOte that all design variables c; , i = 1, 2, ... 6 are almost of the same order, and h(x) and dh(x)/dx are continuous at x 


(4.1) 


We consider a spacecraft with two identical flexible appendages. For simplicity 
made of a single uniform material. 


we assume the appendages are 
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Figure 4 Design of the cross section of appendages 


We begin solving the problem by finding a reasonable number of modes in the reduced order model 
reasonable baseline design with flange equal 4.00 cm. As shown in Table 1, we note that it is appropriate 
flexible modes for a postmaneuver maximum angular deviation to be guaranteed less than 0.05 deg. 


We use a 
to retain three 


Number of modes 
retained in the model 


Max. angle deviation 
post maneuver (deg) 


1.495 0.082 


0.0113 2.6e-3 

Table 1. 


8.3e-4 3.2e-4 1.5e-4 7.1e-5 


Appendage material density, p 

Appendage Material Elasticity, E 
Radius of the rigid central body, R 
Mass of the rigid central body 
Length of one appendage, L 
Maximum torque available, Uo 
Width of the web, b 
Thickness of the web, tl 
Thickness of the flange, t2 
Distributed pay load mass, dm 
Concentrated pay load mass (at x = L), M 
The resource constraint of two aopendages 
The minimal flange depth 
The maximal flange depth 


Spacecraft Data 

3 

1880.00 kg/m 

2.76E11 N/m 2 

12.00 m 

4500.00 kg 

50.00 m 

3.0 E4 N-m 

5.00 cm 
1.75 cm 
0.75 cm 

9.00 kg/m 
None 

450.0 kg 

2.00 cm 

12.00 cm 


Case 1 : Single maneuver case 
Command slew angle, 0 

Thus the exact solution is which is equal to if 

Result : 


90.00 deg 

over the entire feasible design space. 


(a) . 

(b) . 

(c) . 


L l( tf) - 2 V ei = 21.9814 sec, but t f * 3 of this design is 24.6213 sec. 

L A t * ) = tfl = 22 3126 SCC ’ 3nd tf * 3 ° f thiS dcsign is 22 - 41457 sec - switching times between 0 and tr* 3 
of the time-optimal control history arc 1.5547E-8, 0.21945, 0.48124 see (one switching time is almost zero). 

rom <p 2 : tf over domain js 22.3218 sec. Let it be t/ a 
a 5 J 

¥ - L 2 (tf)\ = 9.2E-3. We can accept this design as the solution (as shown in Fig. 5). 


Properties of this optimal design of Case 1 : 
Structural mass of two appendages 
Total pay load mass along the appendages 


379.687 kg 

900.00 kg 
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Total mass of the spacecraft 5779.687 kg 

Total rotational moment 2375330.68 kg-m 2 

Natural frequency <*>i , i = 1, 2, .... 4 0.5642, 1.6942, 4.4738, 8.9745 (rad/sec) 

The max. angle deviation from the uncontrolled modes 0.00908 deg 

Number of switches between 0 and l f /2 of the time-optimal control history : None 

Case 2 : General Multiple Maneuvers 

The set of maneuvers are { 9i } = (9, 15, 30, 45, 60, 90 (deg) }, and assume that they occur at the same frequency. 
Thus the objective function (maneuverability index) is 

«$> = iS'/(o. ) 


i =1 


(4.3) 


* 3 

The solution |d equals M- over the entire feasible design space. 


Result 


(a) . 

(b) . 


L 6 i([L)= 2 V 9i *^o = 13.1753 sec. 


£ — p 1 - 15.0436 sec, and |J. ^ for this design is 15.30617 sec. 


(c) 


(d). 


As : if we let JJ 3 = 15.30617, we have I \X & - L 6 2 ( fljl = 0.26257, as 1.7454 %. 

|4 2 is 14.8580 sec, and |4 ^ for this design is 14.96326 sec. 

As : if we let ^t a = 14.96326, we have I (f 3 - L b 2 ( jx^( = 0.06526, as 0.4392 %. We accept it as the 
solution (as shown in Fig. 5). 

We investigate the exact solution by (P-\ and obtain |I 3 is 14.9455 sec. 

Properties of this optimal design of Case 2 : 

Structural mass of two appendages 425.075 kg 

Total pay load mass along the appendages 900.00 kg 

Total mass of the spacecraft 5825.075 kg 

Total rotational moment 2379168.55 kg-m 2 

Natural frequency w; , i = 1, 2 4 0.8460, 2.0276, 5.5051, 10.6193 (rad/sec) 

The max. angle deviation from the uncontrolled modes 0.02436 deg 

Number of switches between 0 and l f /2 of the time-optimal control history : Three 


V. Conclusion and Future Work 


The problem of combined design of structures and controls for optimal maneuverability of an elastic spacecraft has 
been considered. The main results of the present work are 

i) . The problem formulation is consistent with bang-bang forms of time optimal controls. 

ii) . The performance degradation constraint is considered in the design problem. 

iii) . The optimal design problem is well defined. There always exists a solution. 

iv) . The optimization is done by mathematical programming. 

v) . The gradient of the objective function is computed using the Implicit Function Theorem. 

vi) . Efficient and practical approximate methods have been developed. 

Our experience with various numerical examples leads to the following assertions : 

i). The best structural designs often arc those for which the designs of mass distribution and stiffness distribution have 
very little coupling. 
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ii). The benefit of multiple controls is not apparent, since we can use scalar control to achieve good results. 


Since spacecraft structure is modelled to be linear, with small displacement and inextensible deformation, the 
performace for a realistic system which violates these assumptions is worth investigating. The constraints of structural 
dynamic response, such as maximal displacement and stress, should be considered in the examples as well. Those topics 
are indicated for future study. 



Optimal Design of Appandagaa (Casa 2) 



Position (meter) 


Figure. 5 
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GROUND TEST PROGRAM FOR NEW ATLAS PAYLOAD FAIRINGS 


Introduction 


An extensive ground test program is currently being undertaken by General Dynamics/Space 
Systems Division to verify the design of the metal payload fairings for the new family of Atlas 
launch vehicles, the first of which will be launched this summer. Two new designs, an 1 1-foot and 
a 14-foot diameter version of the payload fairings (see Figure 1), are now available to mission 
planners seeking to accommodate the widest variety of mission requirements. While the 14-foot 
diameter version was developed for the commercial Adas I program, and the 1 1-foot diameter 
fairing was developed for the U.S. Air Force Atlas II vehicle, once both production lines are at full 
capacity, the selection of a fairing will be dictated by the size of the satellite payload. These new 
fairing designs replace the 10-foot diameter, honeycomb fiberglass payload fairings which were 
flown on previous Atlas/ Centaur launch vehicles. The new metal fairings feature a larger payload 
envelope, greater ease of manufacturing and modification, have more consistent quality control 
properties, provide better EMI shielding for the satellite payload, and do so at costs and weights 
comparable to the old fiberglass fairing. 6 


Both the 1 1-foot and 14-foot diameter designs are of aluminum skin, frame, and stringer 
construction and are built at the General Dynamics Services Company plant in Harlingen, Texas. 

1 he mam structural purposes of the payload fairing are to protect the satellite payload during the 
ascent phase and to provide an aerodynamic forward surface for the launch vehicle. After the 
vehicle has cleared the atmosphere, the payload fairing is no longer required, and it is jettisoned 
p lu Sa i V i e ) vei S ht , ar ] c J all ° w f° r the separation of the Centaur upper stage and the spacecraft. 

1 ~ f ° ot ^ nd 1 . 4 : foot designs use a method of separation similar to that originally used for 
the 10-foot fiberglass fairing. At the moment of jettison, which occurs about 3 1/2 minutes after 


STATION 

(INCHES) 


5 0 0 



14-FOOT DIAMETER 11-FOOT DIAMETER 

PAYLOAD FAIRING PAYLOAD FAIRING 

(ATLAS I) (ATLAS II) 

Figure 1. Atlas 11-Foot and 14^oot Diameter Payload Fairings. 



liftoff explosive bolts fire which allow the two 180-degree halves of the fairing to begin 
separation. Spring loaded actuators at the top of the cone section push the halves apart, while the 
aft end of each fairing half begins to rotate on hinges located on the stub adapter. After the fairing 
halves have rotated about 70 degrees, the hinges allow the fairing to safely and completely separate 
from the vehicle. Both fairing halves fall back to Earth, where they land in the Atlantic Ocean. No 
attempt is planned to recover these items. 


Five separate ground tests, one of which has already been completed, were planned to gather the 
necessary data to qualify these new designs for flight. All tests planned will be performed on 
full-scale payload fairing structures (either dedicated test articles or flight articles). Three tests have 
been planned for the 14-foot diameter payload fairing and two tests have been planned for the 
1 1 -foot diameter payload fairing. All tests of the 14-foot diameter payload fairing must be 
complete by the planned June 1990 launch for the Combined Release and Radiation Effects 
Satellite (CRRES) and all tests of the 1 1-foot diameter payload fairing must be complete in time to 
support stress and dynamics analyses which must be performed prior to the January 1991 Initial 
Launch Capability (ILC) date for the Atlas II system. 


T. 14-Foot Diameter Pavload Fairinr Jettison Test 

Description: Of the three ground tests planned for the 14-foot diameter version of the new 
payload fairings, the jettison test is the only one which has been completed as of this writing. This 
test was successfully performed in December 1989 - January 1990 at the Space Power Facility 
(SPF) operated by NASA/Lewis Research Center at the Plum Brook Station near Sandusky, Ohio. 
This site was chosen because it is the largest vacuum chamber in the world, and is the only one in 
which a full jettison of at least one full payload fairing half could be accomplished. The interior of 
the chamber consisted of a metal-walled pressure vessel with a 100-foot diameter circular floor and 
a 120-foot high, domed ceiling. This was surrounded by a thick concrete-walled containment 
building with a profile which betrays the original intent of the building: to house nuclear powered 
satellites during test and checkout. Because the facility had not been used since a 1974 Skylab test, 
a substantial effort was required to reactivate the chamber. Now that the chamber has been proved 
to be operational, several other jettison tests of payload fairings, including one of the giant Titan IV 
fairing, have been planned for the NASA Plum Brook SPF. 

The test article for the jettison test was a dedicated test fairing which was manufactured to the same 
engineering prints and quality standards as a regular flight fairing. This article was the first 
payload fairing completed at the General Dynamics Services Company, Harlingen, Texas assembly 
plant. From its interface with the Centaur upper stage’s stub adapter to the tip of its nose, the new 
fairing stands 39 1/2 feet high. In the SPF, atop its base fixture and stub adapter, the test article 
measured in with an impressive height of 52 1/2 feet. Before each completed payload fairing is 
shipped from its assembly plant, an acceptance test, consisting of a stackmate and a rotation test, 
is conducted to verify that manufacturing tolerances were maintained throughout the entire 
structure. The rotation test consists of splitting the fairing halves (by a total of no more than one 
foot at the top of the cone) using a manually activated screwjack in place of the spring loaded 
actuators. As the fairing halves are gradually rotated on the jettison hinges, clearances at various 
locations along the splitline longerons are recorded to verify all shear pins are disengaging 
smoothly. 

The jettison test was a simulation of the event described above, wherein the two halves of the 
payload fairing separate after the launch vehicle clears the Earth's atmosphere. The test program 
consisted of performing two separate payload fairing jettison events at a simulated altitude of 
85,000 feet (chamber pressure = 17 torr). After the inner metal door and the outer concrete 
chamber doors were closed, an approximately five hour pumpdown was performed, and after a 
short countdown a switch was flipped, immediately supplying power to the explosive bolts, 
initiating the jettison event. For each event, one half of the fairing (capped half) was fully 
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jettisoned, while the other half rotated only fifteen degrees before impacting the catch net (see 

Sr 2 u- B i° th halves were ? lowed from initial net im pact to a full stop by hydraulically actuated 
brakes which were an integral part of the catch net system. The force used to initiate the jettison 
event was varied by using only one spring-loaded jettison actuator for the first event and J the full 
complement of two actuators for the second jettison event. This provided analysTwkh two data 
&r nSt WhlCh t0 compare 1116 analytically predicted behavior of the test fairing After the 

Xrsr oniy a few days werc » » » “7 

primar y Purpose of the jettison test was to demonstrate that the analytical 
^ ST kT N c . omput . er model being used by structural dynamicists and stress analysts to predict 
airing behavior dun ng vehicle flight environments is able to accurately predict the behavior of the 

Und r r the «? n test conditions. Pre-test predictions of all test dSJtoSShoS 
gid body motion, fainng half-breathing modes) were made using the analytical computer model ’ 

A u \? ese anaIy,ical prcdictions wil1 indi “ le if any SSL ™ e • 

necessary to the model Another purpose was to simply demonstrate that the fairing jettison 
hardware (actuators hinges, explosive bolts, shear pins, harness disconnects e™) Vuncrioned 

bet^Sn Pf!iL A h n important element of proper jettison functioning is the mechanical clearance 
een the fainng hardware as it rotates on the hinges and other critical items on the Centaur 
upper stage and the satellite payload. Because the fairing half "breathing” or "pinching" mode 
effectively reduces the static clearance, the intent was to design a fairing whW was as stiTf a, 

, The ll 9 ms whlch represented potential rotation interferences were simulated during this 

b ankets for intenor noise reduction, were installed inside the fairing test micle to verify that the 
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Figure 2. Payload Fairing Jettison Test Catch Net System. 


ORIGINAL PAGE IS 

OF POOR QUALITY 



bending associated with the jettison event would not damage either the mylar components or their 
fastening hardware. It was felt that this verification was necessary to eliminate the possibility that 
loose mylar could damage the encapsulated spacecraft during the jettison event. 

Data Acquisitions; Test instrumentation for the jettison test included 31 channels of 
low-frequency accelerometer data, intended to measure the lowest vibration modes of the fairing 
halves and larger scale, rigid-body motions, 21 channels of high-frequency accelerometer data, 
which measured the shock environment (generated by the 28 explosive bolts) at various sensitive 
locations on the vehicle structure, and 44 channels of strain gage data, which measured the 
structural loads in the jettison hinges on which the fairing halves rotate. An analog FM-based data 
acquistion system was used to record all test data. Events were recorded by fourteen high-speed 
motion picture cameras and three video cameras. The video cameras were installed mainly to 
present to the test conductor a real-time image of what was occurring inside the chamber, realizing 
that there are no windows through which to look. Exactly half of the fourteen high-speed film 
cameras were focused on clearly visible targets mounted on fairing hardware. This film was later 
run through a film motion analyzer to produce deflection, velocity, and acceleration data 
Clearance indicators, made of thin solder wire, were used to detect any infringement of the payload 
fairing structure during its jettison rotation into sensitive areas surrounding avionics packages or 
into the satellite payload envelope. 


TT. 14-Foot Diameter Pavload Fairing S tructural Test 

Description: This test is currently under way at the Sycamore Canyon test site operated by 
General Dynamics just north of San Diego, California. During this structural test, the fairing will 
undergo a series of static test conditions which will determine if the fairing structure will yield at 
predicted design limit loads (based on loads experienced during the transonic condition) or will fail 
at design ultimate loads (125% of design limit loads). Because the test article is to be exposed to 
ultimate structural loads, the dedicated test fairing which was used for the jettison test described 
above, and which will never fly on an Atlas vehicle, will be used for this test. Four different test 
configurations (see Figure 3) will be used in order to completely test all of the major structural 
elements of the payload fairing: upper and lower nose cone (crush pressure and side load), nose 
dome (crush pressure only), cylinder and boattail (burst pressure and vent fin loads), and the all-up 
system level configuration (bending moments, shear loads, and axial loads). The last two 
configurations will be performed in a new test tower constructed last year specifically for the 
purpose of performing payload fairing testing. Included as part of the all-up system level testing 
will be test conditions which will reveal stiffness data on the payload fairing structure. It is 
planned that 27 separate test conditions will be required to fully accomplish the objectives of this 
test. The following paragraphs describe the major test configurations and the purpose for each of 
these: 

Nose Cone Test Conditions: For the nose cone tests (see Figure 3a), the cone section (21 
feet high) will be removed from the fairing cylindrical section and mounted on an airtight base 
fixture. A negative pressure differential will be established across the nose cone skin which will 
simulate worst-case crush pressures experienced by this structure during the ascent phase of flight. 
While the fairing structure is vented at the bottom of the cylinder section, pressure differentials are 
still experienced during flight at different fairing stations due to the varying aerodynamic pressure 
profiles. Pressures in the cone region during vehicle ascent are of the crush variety due to the 
aerodynamic nature and purpose of this structure. A shear load will also be introduced at the top of 
the cone during these test conditions in order to observe and characterize post-buckling behavior 
and load carrying capability of the monocoque (no external stringers) cone structure. It should be 
mentioned that buckling of the structure is expected at high loading conditions and is not to be 
considered a failure of the structure. There will be a total of four test conditions devoted to nose 
cone testing. 
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Figure 3a. 

Nose Cone Crush Tests. 
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Figure 3b. 

Nose Dome Crush Tests. 


internal pressure 
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Cylinder Burst Tests. 


NOSECOf€ 
AXIAL LOAD 



Figure 3d. 

All-Up System Tests. 


Figure 3. Payload Fairing Structural Test Configurations. 
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Nncp Dome Tp^ Conditions: The nose dome crush pressure test conditions are fairly straight- 
forward? involving only the creation of a negative pressure differential across the nose dome 
structure, simulating the crush pressures seen during the ascent phase of the vehicle s flight (see 
Figure 3b). No shear, bending, or other structural loads will be imparted into the nose dome 
during these tests. There are only two nose dome test conditions. 

C ylinder and Boattail Test Conditions; A unique test fixture was required for the cylinder 
and boattail burst pressure test conditions (see Figure 3c). Pressurizing the entire volume of the 
cylinder using a flat disk to seal the top of the cylinder (about 22,100 square inches) was not an 
option because this would have imparted large axial tension loads into the cylinder skin. Instead, a 
cylinder of a slightly smaller radius will be inserted inside the fairing cylinder, and the small 
annulus between the two cylinders will be pressurized. The primary goal in these test conditions is 
to observe the behavior of the explosive bolts and the split line longerons. Gapping or t e 
longerons in the areas between bolts will be measured to characterize this behavior. No shear 
loads will be imparted to the cylinder or boattail structure during this test, but there will be small 
loads input into the vent Fin structure to observe how flight loads from the vent fin are distributed 
into the cylinder skin and the backing frames and stringers. There will be a total of four cylinder 
and boattail test conditions. 

All-IJn System Level Test Conditions*. There are 17 test conditions to be conducted in the 
fully assembled, all-up test configuration (see Figure 3d), including five stiffness test conditions. 

In this configuration, die test article will be assembled in exactly the same way that the fainng will 
sit on top of the Centaur upper stage during the boost phase of flight. Axial loads will be imparted 
both at the top of the nose cone and at the boattail/ cylinder interface, and side loads will be able to 
be input at two different stations on the nose cone. The 12 non-stiffness test conditions are 
constituted by combining three test configurations in four variations: using design limit loads, 
using design ultimate loads, loading parallel to the split line, and loading perpendicular to the split 
line. The three configurations are side loading at the upper fixture on the cone (which tests the 
upper portion of the cylinder), side loading at the lower fixture on the cone combined with . 
maximum axial compression loads, and side loading at the lower fixture on the cone with minimum 
axial compression loads (both of which test the lower portion of the cylinder, the boattail, and the 
stub adapter). The maximum axial load condition imparts worst-case compression loads into the 
skin, longerons, and explosive bolts on the compression side of the fairing, and the minimum axial 
load' condition imparts worst-case tension loads into the skin and longerons on the tension side of 
the fairing. 


Test Instrumentation: About 300 strain gage channels and about 45 deflection transducers will 
be present on the test article. The strain gage locations are distributed evenly on the various 
components of the test article structure, but typically only half of the gages will be read for a given 
test condition. The strain gages are intended to provide an indication of how loads are distributed 
throughout the test article structure. Deflection transducers are present mainly to provide data on 
the stiffness of the test article and to monitor deflections of the test article to ensure safe test 
operations. Load cells on test fixtures and pressure transducers in the test setup will also be 
monitored to obtain an accurate picture of loads and pressure differentials being input into the test 
article. A set of digital data loggers will be used for data acquisition system, and data will be 
delivered to the stress department in personal computer compatible spreadsheet formats. 
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IU. 14-Foot Diameter P avload Fairing Acoustic Test 

Description: The acoustic test of the payload fairing will be conducted during May 1990 in 
General Dynamics/Space Systems Division's Acoustic/Thermal Test Facility (ATTF) at the Kearny 
Mesa plant in San Diego, California. The ATTF is a dual chamber containing separate acoustic and 
thermal test chambers. Both chambers are maintained as 100,000 class clean facilities. Its acoustic 
chamber is one of the largest acoustic test facilities with a floor measuring 33 x 40 feet and a ceiling 
height of 50 feet (65,000 cubic feet). The chamber is fully reverberant, having a 25-Hz horn, a 50- 
Hz horn, and two 100-Hz cutoff horns, all mounted on the north wall of the chamber (see Figure 
4). The frequency range of the chamber is 25 - 10,000 Hz, and the rated overall sound pressure 
level is 154 dB. Chamber environments and data acquistion are controlled from the control room 
on a mezzanine above the chambers. The primary purpose of the ATTF is the environmental 
testing of large space structures. Currently, the major emphasis of the chamber is on the rigorous 
checkouts required prior to delivery of the Centaur upper stage of the USAF Titan IV launch 
vehicle. This checkout includes a full exercise of Centaur avionics control systems while the 
vehicle L02 tank is filled with liquid nitrogen and subjected to launch acoustic levels. The vehicle 
with empty propellant tanks was later subjected to thermal cycling ranging from -40 F to 185 F. 

During the acoustic test of the Atlas 14-foot diameter payload fairing, a fully assembled fairing will 
be subjected to acoustic levels representative of both the launch and Max Q environments. External 
sound pressure levels for both conditions will approach the rated capacity for the chamber. Empty 
chamber calibrations will be performed prior to arrival of the test article to better characterize the 
obtainable sound pressure levels. Various re-configurations of the test article will be performed in 
order to characterize an 1.) empty, generic fairing, 2.) a generic fairing with a satellite payload, 3.) 
an empty fairing with the acoustic blanket installation, and 4.) a fairing with the acoustic blanket 
installation and a satellite payload. Other design features which will be tested are the noise 
mufflers which are placed over the vent holes in the fairing cylinder section. The flapper doors in 
the muffler structure (designed to prevent payload contamination and noise intrusion but still allow 
for venting of internal burst pressures) will be alternately opened and closed to determine their 
effect on internal noise levels. The noise reduction properties of the fairing structure and the 
acoustic blanket installation will be fully tested by the completion of the ten test conditions planned 
(8 at launch levels and 2 at Max Q levels). In order to fully characterize the fairing acoustics, 
sound decay measurements will be taken inside the fairing structure prior to testing, both with the 
acoustic blankets installed and with them removed, to determine the reverberant component of the 
measurements which will be taken internal to the fairing during the actual testing. 

Also of interest to the structural dynamicists are the vibrations induced in the fairing structure by 
this acoustic energy. To investigate this phenomena, accelerometers will be placed on mass 
simulated avionics packages mounted on the forward end of the Centaur upper stage and on areas 
of payload fairing skin. By placing microphones very near the avionics packages and the payload 
lairing skin on which accelerometers are placed, data are gathered which will allow transmissibility 
studies to be performed. The vibration environments, while significant for the purposes of this 
test, are expected to be mild enough to allow the use of flight hardware for this test. The second 
flight payload fairing and the third 14-foot diameter article to come off the line at Harlingen will be 
used for this test in order to help acheive some schedule compression (the dedicated test article will 
be committed to the structural test at this time), and because the standard cork and paint installation 
which was incompatible with the goals of the structural test, is a requirement for the acoustic test. 

It was felt that the cork and paint on the nose cone section, bonded on as an ablative for thermal 
control, would have a significant impact on the acoustic transmissibility of the nose cone skin. 

Pflttl AyqpigifiQn; Acoustic levels internal to the payload fairing structure will be fully 
characterized through the use of about 20 microphones positioned both in and around the test 
article. Several control microphones will also be used to monitor chamber conditions in real time. 
Vibrational levels associated with the acoustic energy will also be recorded by about 10 triaxial 
accelerometer placements on sensitive areas of vehicle structure. A state-of-the-art data acquisition 
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Figure 4. General Dynamics Acoustic Chamber, San Diego, California, 




system in the ATTF control room will record all test data and will allow for a wide variety of data 
presentation formats; 1/3 Octave Band Center Frequency Plots shall be produced for all 
microphones. Power Spectral Densities (PSDs) will be produced for all data microphones and 
Accelerometer Spectral Densities (ASDs) of all accelerometer data will also be produced as will 
ttansmissibihty plots showing the relationship between microphone and accelerometer data. 
Because this test is currently scheduled for completion less than one month prior to launch of 

AC-69, the first commencial Atlas flight, prompt data reduction and presentation will be of utmost 
importance. 


Descr i ption ; This first test of the smaller 1 1 -foot diameter fairing, developed for the USAF 
Atlas II, will be conducted in two major segments: a "free-free" test condition in which only one 
half of the fairing will be tested, and a "fixed base" test condition which will be perfonned on a 
fully assembled and erected fairing. The free-free conditions will be performed at the General 
Dynamics operated U. S. Air Force Plant 19 near downtown San DiegQ and the fixed base test 
conditions will be Performed at the Sycamore Canyon test facility. The fixed base testing will 
immediately follow the 14-foot diameter payload fairing structural test and be conducted^ the 

T™: ‘ e n S “ n ° Wer - The test a J lcle .y i11 be a dedicated 1 1-foot diameter test article, constructed to the 
same engineering prints and quality control criteria as a flight article. This will be the first 1 1 -foot 

™lti?nH r / ainn8 c ° m P Ieted . by ^ Harlingen assembly plant and will be the fourth aluminum 
P ay l° d f g offth e Production line. Because of suspicions that the thermal control cork coating 
H ?? ne W0U d ; ntrodu f sl S mflcant modal damping, it was determined that it should be g 
llfnmH° n h f test T 1C l e f in ° rder t0 accuratel y reflect the properties of the flight article The 
11 -foot diameter payload fainng is somewhat shorter than the 14-foot diameter version measuring 

the nsrs m ght ’ 5 1/2 / ee J le ?5 V han its larger brother - WhiIe * is enough ^accomSate g 
!he s P acecraft > the 1 T-foot diameter payload fairing is some 1,500 pounds lighter than 

e 14-foot diameter payload fainng, weighing in at just over 3,000 pounds. This allows for a 
dramatic increase in payload-weight-to-orbit. The 1 1-foot fairing is even slightly lighter than the 
comparably sized 10-foot diameter fiberglass fairing, showing the inherent efficfency of the 
aluminum skin stringer, and frame construction. In addition to supplying data on fairing 

be the prim ^ — ° f 

experience in t£e ^ C ° n,raCt age " Cy ' SDRC> WWch haS a ^ deal of 

• J!!* 01 ™ a P re ' test analyses, which will predict mode shapes and frequencies 
determine instrumentation quantities and placements 

• install the instrumentation and shakers on the test article 

• perform all test operations 

• perform all data acquisition 

• perform a post-test analysis intended to refine the fairing analytical model 

•prepare a test report which presents all test data yi 

a -gyFrte Tfst Emimnmla; The free-free test condition, in which a sinele fa inn,, half will h. 
exctted ts being petformed todetenoine the vibrational modes of a fa Mn E 
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fairing will be suspended from the ceiling in a concave down (inverted canoe, see Figure 5) 
orientation and will be supported at three locations by bungee cords. Exciters will be mounted on the 
floor and will be attached to the split line longerons. Multiple excnauo J will be the 

primary method of excitation, providing information on all modes below 50 Hz, but for the major 
modesof interest, sine-sweep excitadon will be used to obtain more specific information Ojnean y, 
orthogonality damping, etc.)* Orthogonality requirements for the free-free condition are that all 
off-diagonal terms of the modal mass matrix be less than 0. 10. Up to 12 total retakes of contaminated 
target modes will be permitted to satisfy these requirements. 

Fixe d Rasp Test Requirements: In order to determine the vibrational modes of the fully 
assembled fairing as it sits on the Centaur upper stage for the first 3 1/2 minutes of flight, the fbced 
base modal configuration will be performed at Sycamore Canyon. The modes of primary interest are 
the first three bending modes in the vehicle pitch and yaw axes, the first two axial modes and the first 
torsional mode. If any of these modes is above 64 Hz, a shell mode pair will be substituted. It is 
anticipated that about 250 low-frequency accelerometers will be required to monitor vibration during 
this test condition. As in the free-free condition, multiple input random excitation will be the primary 
source, characterizing all modes below 50 Hz, and sine-sweep will be used to isolate the .speeded 
target modes. Orthogonality requirements are, as in the free-free condition, that all off-diagonal terms 
of the modal mass matrix be less than 0.10. A total of 18 retakes will be allowed in order to satisfy 
these requirements for all contaminated target modes. 




view A - A 


Figure 5. 1 1-Foot Diameter Payload Fairing Modal Test Free-Free Condition. 
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■Y> 11-Fwt Piamgter Pavload Fairin g Structural Tcct 

D f S CnptlQn; The 1 1-foot diameter payload fairing structural test will be the last fairing test to be 
conducted in the Sycamore Canyon test tower. There will be no special setup requhS fonhis teVt as 

for t?!E° n and t th ! te l ai i! 1Cle f ° r . this te , st exactl y the same as the fixed base configuration 
r f hn,^^TnM^ht U H f ey f r5 1Ch ?• W1 imme ^ately follow. The goals of this test are very similar to 
ln?viS abl ‘ she ^ for th . e lf-foot diameter payload fairing structural test: 1.) the fairing structure shall 
°J ^ eld under design limit loads and 2.) the fainng structure shall not fail under design ultimate loads 
iinH v of d ® sl ? n h ™ u) - As with the 14-foot diameter structural test, there will be a series of static test 

sriffop^flh 518116 ? S Pf C ; f,Cally l ° sa0sfy the above requirements and additionally to determine the 
stiffness of the payload fainng structure, serving as a verification of the stiffness data obtained during 
the previously mentioned modal survey test. This structural test, however will not incorporate anv of 
die component level tests which the 14-foot structural test program did, mainly because Tome of the 
structural components tested at that time (upper and lower nose cone, nose dome, and stubTLpter) 

^e common between the 1 1-foot and 14-foot designs and need not be demonstrated again The only 
guration for this test is the fully assembled fairing. No burst or crush pressures will be required 
during any of the test conditions to be performed in this test program. ^ 

X sst Conditidnfl; Like the all-up system level tests conducted during the 14-foot diameter 

^ C h a J, teSl ? ro S ram ’ , oadin .£ of the 1 1 ‘ foot diameter test article will be accomplished through a 
combination of axial and side loading. There will be axial loading available through a fixture at the 
top of the nose cone and side loading available at three locations: the top of the nose cone near the 

mrrempmal 11 ?? and near the mid dle of the fairing cylinder section (see Figure 6). By inputing 
incremental side loads at three stations, a very good approximation of flight shear and bending P § 

moment profilescan be obtained, avoiding any seriously overloaded structure. A total of eleven test 
conditions will be required to complete the structural test program. There will be four standard 
conditions using maximum axial compression loads combined with design limit and ultimate shear 


THREE SIDE 
LOAD FIXTURES 



Figure 6. 1 1-Foot Diameter Payload Fairing Structural Test Loading Fixtures. 


loads/bending moments conducted both parallel and perpendicular to the fainng split line. .Ttej wiU 
also be two test conditions which use mimimum axial compression loads in conjunction with side 
S (5S and ultimate) parallel to the fairing split line. Maximum axtal C™P = relays 
to the Max Q flight condition and the minimum axial compression case corresponds to loading o 
Atlas vehicle at *e transonic flight condition. Two test conditions (design limit and ultimate loads) 
tin tCim that Ae axis of compression be off the standard axes which ate paral el/perpendK inter I to 
the spht line. This is required because of an air conditioning duct door which is located 30 decree 
off the split line. Because of the location of this door at the very top of the fainng cylinder, this 
smicture P cannot be designed to a factor of safety of 2.0, and as such must be tested to ultimate loading 
conditions per the USAF contract. An axis of compression which is off the standard axes squires 
that a resultant load be reacted, dictating that twice as many load cells be used to impart the test loads. 
The final three test conditions are the stiffness tests, one each in the standard side load axes (parallel 
and perpendicular to the split line) using pure side load, and one pure axial load condition. 

Data Acauisition: Test instrumentation for the 1 1-foot diameter payload fairing structural test 
includes about 80 strain gages and about 35 deflection transducers. Strain gages will be placed in 
circumferential patterns at the bottom of the lower cone and at the mid-point of the small transition 
cone between the cylinder and the split barrel. Several explosive bolts shall also be instrumented with 
strain gages. Deflection transducers will be used mainly to collect stiffness (deflection versus load) data 
and to monitor test safety during high loading conditions. Load cells connected to the hydraulic load 
cylinders will also be monitored (maximum of seven hydraulic load cylinders operating during 
off-axis air conditioning door test conditions) to verify the loads input into the test article. The data 
acquisition system used on this test will be identical to that used during the 14-foot diameter payload 
fairing structural test program. A digital data logger will record data on the 3 1/2 -inch floppy disks in 
standard personal computer spreadsheet format. 


Conclusions 

To establish the competitiveness of the revitalized family of Atlas launch vehicles (I, II, IIA, and 
IIAS) a new series of payload fairings, an 1 1-foot and a 14-foot diameter version, were designed to 
accommodate the widest possible variety of satellites. Because these aluminum fairings are new 
designs, the plant at which they are produced is new, and launch customers are very 
anxious to fly their payloads, an ambitious and efficient test program is essential. Five major tests 
have been planned for completion within the span of one calendar year. One of these has been 
completed, with every indication that it was a success, one is currently under way, and two more are 
scheduled to start in the month of April. Through effective use of test assets, facilities, and 
personnel, all testing will be completed, allowing the fairing design to be completely characterized and 
then qualified through analysis prior to first launch of each of the fairings. 
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DIRECT USE OF LINEAR TIME-DOMAIN AERODYNAMICS IN 
AEROSERVOELASTIC ANALYSIS: AERODYNAMIC MODEL 


J. A. Woods and M. G. Gilbert 
Aeroservoelasticity Branch 
NASA Langley Research Center 
Hampton, Virginia 


PRECEDING 
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The work presented here is the first part of a continuing effort to expanding existing 
capabilities in aeroelasticity by developing the methodology which is necessary to utilize 
unsteady time-domain aerodynamics directly in aeroservoelastic design and analysis. 

The ultimate objective of this study is to define a fully integrated state-space model 
of an aeroelastic vehicle's aerodynamics, structure and controls which may be used to 
efficiendy determine the vehicle's aeroservoelastic stability. 

In this presentation, the current status of developing a state-space model for linear 
or near-linear time-domain indicial aerodynamic forces is presented. 


MOTIVATION: 

TO EXPAND EXISTING AEROSERVOELASTIC DESIGN AND 
ANALYSIS CAPABILITIES TO INCLUDE THE USE OF 
UNSTEADY TIME-DOMAIN AERODYNAMICS 

LONG-TERM OBJECTIVE: 

DEVELOP METHODOLOGY TO UTILIZE LINEAR AND NEAR- 
LINEAR TIME-DOMAIN AERODYNAMICS IN THE SUPERSONIC 
AND SUBSONIC REGIMES DIRECTLY IN AEROSERVOELASTIC 
DESIGN AND ANALYSIS. 

IMMEDIATE OBJECTIVE: 

DEVELOP A TIME-DOMAIN STATE-SPACE MODEL OF TIME- 
DOMAIN AERODYNAMIC INDICIAL FORCES. 
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THE INTEGRATED AEROELASTIC MODEL 

To understand the importance of this research, it is necessary to consider that 
several codes [1,2] have been developed in recent years which compute time-domain 
unsteady aerodynamics, however, the techniques needed to utilize the aerodynamics in 

aeroservoelastic design have not been fully developed. 

One of the only methods devised to date to evaluate the aeroelastic stability of 
aerospace vehicles in the time-domain has been a general method capable of handling the 
nonlinear system [3]. This method is expensive as it involves the computation of the 
aeroelastic system time response which requires solution of the nonlinear small disturbance 
aerodynamic equations. Further, a frequency decomposition of the response is necessary 
to evaluate the stability of component modes. The response must be recomputed at several 
dynamic pressures until a neutrally stable mode is encountered. Other available methods 
model the aerodynamics directly in the frequency domain. 

For linear and near linear systems in supersonic and subsonic flow, however, the 
vehicle stability may be evaluated without computing the aeroelastic system forced response 
or transforming forces to the frequency domain. This is accomplished by representing the 
time-dependent aerodynamic forces in state-space form coupled with a commonly used 
state-space representation of the structure. Stability is determined by the eigenvalues of the 
coupled system matrix. 

The focus of this presentation is, again, on the formulation of the aerodynamic 
portion of the integrated model. 



Figure 1. Schematic block diagram indicating integration of the aerodynamic 
model with the structural model. 
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FORMULATION OF AERODYNAMIC MODEL 


The aerodynamic model is derived as the Laplace transform of a commonly used 

frequency domain approximation modified from ref. 4. It is transformed directly into state- 
space form. 


MODIFIED FREQUENCY DOMAIN APPROXIMATION 
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APPROXIMATION METHOD 


The approximation method involves a least squares approximation to the actual 
aerodynamic force to determine the scalars A 0 , Aj and Bj. The fit is constrained at t=0 to 

fit exactly and at large times to equal the asymptotic value of the generalized force. As in 

the frequency-domain rational function type approximations, aerodynamic poles, P j , are 

initially specified. 

The aerodynamic forces currently being approximated are the rigid-body forces 
acting on a NACA0064 airfoil and are due to Dowell [5]. 


APPROXIMATING FUNCTION 
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APPROXIMATION METHOD (CONTI) 


A system identification technique frequently used in control system analysis is 
applied to regenerate the generalized aerodynamic force. Specifically, impulse and step 
responses of the aerodynamic model are generated using a discrete-time state-transition 
method. The sum of these responses is the aerodynamic approximation, Q(t)=Q(t) , based 
on previously determined coefficients and the specified aerodynamic poles. 

Due to the discontinuity at t=0 in the impulse input, an assumption is made that at 
t=0+, initial conditions are real valued. At t=0-, initial conditions are zero. This 
assumption can be shown mathematically. 


STATE TRANSITION EQUATIONS 

w(t + i) = ow(t) + ru(t) 

Q(t) = Cw(t) + Du(t) 

WHERE 

*(t) = e fA]T and r = /Je [A]t Bdx 

ASSUMING 

w(<)-) = 0 w(0 + ) = Bu(0) 
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APPROXIMATION METHOD (CONTI) 


Improvements to the aerodynamic approximation are made by updating the 
aerodynamic poles, Pi, followed by another least squares approximation to recompute the 
coefficients. To update the poles, the method used by Peterson and Crawley [7] to 
approximate unsteady aerodynamics in the frequency domain is implemented in the time 
domain. A norm square-error cost function is defined. In this case, the square of the 
difference between the actual aerodynamic force and the approximation is used.. The 
incremental change in aerodynamic poles is solved for by inverting the Hessian, 

d 2 J / dPjdpk, in a single term Taylor series expansion of 3J / 3p. The incremental change 
in Pj is multiplied by a scale factor, a, and added to the current aerodynamic poles. The 
scale factor, a, is computed using quadratic interpolation [8] to insure that the cost is 
approaching a local extrema. 

The new aerodynamic poles are limited. If a given pole is greater than -0.01, it is 
set equal to that value until the next parameter update. To prevent a pole from going to -<*> 
and ill-conditioning the system matrix later on, the pole is limited to a value which would 
produce no more than a 99.5% decrease in magnitude of the exponential over a given time 

step. 

The two step procedure of computing system coefficients and updating 
aerodynamic states is repeated until the cost function has been minimized. 


SQUARE ERROR COST FUNCTION 

J(p)=[Q(t)-Q(t,p)] T [Q(t)-Q(t,p)] 
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AERODYNAMIC POLE UPDATE 
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PRELIMINARY RESULTS 


Table I briefly describes some of the progress which has been made up to this time. 
Four sets of initially specified aerodynamic poles, associated coefficients and the initial cost 
are indicated as well as the minimum cost quantities. The first set of poles is a subset of the 
poles which were used by Dowell to generate the aerodynamic forces. Dowell's zero pole 
was not included for stability reasons and because the A\ term serves the same purpose of 

providing a constant term at t=0. The other sets of poles represent "random" selections 
between a small negative number and -1.0, -2.0 and -3.0. 

A minimum cost was obtained for each of these sets of poles. The poles close to 
those of the generating function produced the lowest cost. Minimum cost increases from 
there as the range of initial poles widens. It is noted that finding a minimum isn’t always 
guaranteed. For some sets of initial poles, the least-squares fit doesn't converge or the 
program determines a local maxima instead of a local minima. 

One of the immediate observations which can be made from Table I is that the 
aerodynamic poles tend to decrease in magnitude as the cost is minimized. The same trend 
occurs as other rigid body forces are being approximated. The implication is that the fit 
improves at large times and degrades at small times. In terms of reduced frequencies, this 
means that the high frequency components of the curve are not being fit well. Thus, a 
weighted least-squares fit and a weighted square error function will be considered to 
improve the approximation at small times. 

Finally, an assumption made in the quadratic interpolation subroutine which 
computes parameter step size is that when the square-error cost is computed for the "step- 
ahead coefficients remain constant . From Table I, this appears to be a valid assumption, 
as over the entire range of parameter updates, coefficients have remained fairly unchanged. 
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PRELIMINARY RESULTS 


LIFT DUE TO PLUNGE 
6 POLE APPROXIMATIONS 


Initial Data 


Poles* Coefficients A 


Minimum Cost Results 


Poles Coefficients A 


- 0.1 

- 0.433 

- 0.3 

- 0.51 

- 0.8 

1.646 

- 1.2 

- 2.841 

- 1.75 

2.19 

- 3.5 

- 0.574 

- 0.2 

- 3.041 

- 0.4 

43.06 

- 0.5 

- 124.97 

- 0.6 

121.388 

- 0.8 

- 51.012 

- 1.0 

14.162 

- 0.334 

- 6.036 

- 0.668 

46.844 

- 1.0 

- 178.225 

- 1.334 

333.826 

- 1.668 

- 299.231 

- 2.0 

102.678 

- 0.5 

- 9.818 

- 1.0 

89.054 

- 1.5 

- 365.568 

- 2.0 

723.748 

- 2.5 

- 678.553 

- 3.0 

241.917 


0.999 0.034 0.0170 


0.999 - 0.074 0.0047 


0.999 - 0.343 0.1995 


0.999 - 1.266 0.8607 


- 0.091 

- 1.24 

- 0.081 

0.805 

- 0.148 

- 0.178 

- 0.520 

- 0.116 

- 1.255 

0.402 

- 2.082 

- 0.123 

- 0.185 

- 3.202 

- 0.33 

42.016 

- 0.401 

- 114.617 

- 0.470 

103.931 

- 0.613 

- 37.816 

- 0.775 

9.256 

- 0.243 

- 5.681 

- 0.424 

42.426 

- 0.590 

- 152.419 

- 0.747 

263.921 

- 0.898 

- 215.887 

- 1.049 

67.253 

- 0.35 

- 6.406 

- 0.683 

42.793 

- 1.03 

- 137.988 

- 1.409 

222.239 

- 1.837 

- 174.33 

- 2.343 

53.628 


0.999 - 0.036 0.0012 


0.999 - 0.056 0.0020 


0.999 - 0.0997 0.0166 


0.999 


- 0.421 0.2469 


Poles indicated are B; — . 

1 V 


Table I. Summary of some aerodynamic poles, coefficients and cost functions 
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PRELIMINARY RESULTS (CONTI) 


Figure 2 illustrates some of the approximations to the rigid-body forces acting on a 
NACA0064 airfoil which are currendy being obtained. The figure includes a pair of 
figures for each of four rigid-body aerodynamic forces. The lower figure in each pair 
contains a comparison between the aerodynamic data and an approximation made by using 
the initially specified aerodynamic poles. The norm square-error cost is indicated. The 
upper figure in each pair indicates the improved approximation after the minimization 
technique has been applied. Again, the minimum norm square-error cost is indicated. In 
all cases, aerodynamic data has been normalized with the largest absolute magnitude of 
force. 

As can be seen, the technique does improve the approximation noticeably. In three 
of the four cases, the cost has been reduced by about 90%. In the case of "Moment Due to 
Pitch , the cost was observed to remain high even after minimization. This emphasizes the 
fact that the current method finds only the first extrema in cost. This extrema may be only a 
local extrema and not a global one. 
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PRELIMINARY RESULTS 

NACA0064 RIGID BODY AERODYNAMIC FORCES 
P 0 =[-.334 -.668 -1.0 -1.334 -1.668 -2.0] 


Lift Due to Plunge 


Jmin =0.0166 


Lift Due to Pitch 


Jmin = 0.0185 


Jo = 0.1995 


Moment Due to Plunge 


Jmin = 0.1539 


Moment Due to Pitch 


Jmin = 5.5157 


V 


Jo =1.1 159 


Figure 2. Approximations to NACA0064 airfoil rigid body forces using initial aerodynamic 
poles and aerodynamic poles computed for minimum cost. 













FURTHER DEVELOPMENT 


Other methods will be considered to determine minimum cost. The method 
currently used is effective, but needs modification. 

In an effort to improve the fit for small times, a weighted least squares fit will be 
implemented to determine the coefficients. A weighted square error cost function will also 
be considered. 

Sometimes the program converges to a local maximum instead of minimum. Thus, 
means of forcing the program to converge on a minimum will be implemented. 


FURTHER DEVELOPMENT 

• IMPROVE PROCEDURE TO IDENTIFY MINIMUM COST 

• INVESTIGATE WEIGHTED LEAST SQUARES APPROXIMATION TO 
DETERMINE COEFFICIENTS 

• INVESTIGATE A WEIGHTED SQUARE ERROR COST FUNCTION 

• INVESTIGATE METHODS OF CHANGING THE SEARCH DIRECTION IF A 
MAXIMUM IS BEING APPROACHED INSTEAD OF A MINIMUM 
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FURTHER APPLICATIONS 


To further evaluate this technique, aerodynamic data generated for a real aircraft by 
a time-domain aerodynamic code in the subsonic and supersonic flight regimes will be 
modeled. Both rigid-body and flexible modes will be considered. 

Finally, to fulfill the whole purpose of developing this model, methodology will 
need to be developed to integrate the aerodynamic model effectively with a structural 
model. Later, control systems will be integrated into the scheme. Using the integrated 
models, system stability will be evaluated. 


FURTHER APPLICATIONS 

• APPLY TECHNIQUE TO FLEXIBLE AND RIGID BODY GENERALIZED 
AERODYNAMIC FORCES ACTING ON A REAL AIRCRAFT 

• DEVELOP METHODOLOGY FOR INTEGRATING MODEL WITH DISCRETE- 
TIME STRUCTURAL MODEL AND PERFORMING STABILITY ANALYSIS FOR 
ARBITRARY MOTION 
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Motivation 


Time-correlated gust loads are time histones of two or more load quantities due to the same disturbance time 
history. Time correlation provides knowledge of the value (magnitude and sign) of one load when another is 
maximum. At least two analysis methods have been identified (references 1 and 2) that are capable of 
computing maximized time-coirelated gust loads for linear aircraft. Both methods solve for the unit-energy 
gust profile (gust velocity as a function of time) that produces the maximum load at a given location on a 
linear airplane Time-correlated gust loads are obtained by re-applying this gust profile to the airplane and 
computmg multiple simultaneous load responses. Such time histories are physically realizable and mav be 
applied to aircraft structures. y 

Within the past several years there has been much interest in obtaining a practical analysis method which is 
capable of solving the analogous problem for nonlinear aircraft. Such an analysis method has been the focus 
of an international committee of gust loads specialists formed by the U. S. Federal Aviation Administration 
?,.^ as the topic ofa panel discussion at the Gust and Buffet Loads session at the 1989 SDM Conference in 
obile, Alabama. The kinds of nonlinearities common on modem transport aircraft are indicated in figure 1 

The Statical Discrete Gust method of reference 1 is capable of being, but so far has not been, applied to 
pr^edureTsTsfential “ referenCe lf t0 make the method practical for nonlinear applications, a search 


The method of reference 2 is based on Matched Filter Theory and, in its current form, is applicable to linear 
systems only The purpose of the current paper is to present the status of an attempt to extend the matched 
, er approach in reference 2 to nonlinear systems. The extension uses Matched Filter Theory as a starting 
point and then employs a constrained optimization algorithm to attack the nonlinear problem 


• Time-correlated gust loads are generated to obtain physically 
realizable design loads for the analysis and design of aircraft 
structures 

• Active control systems of modern aircraft contain significant 
nonlinearities: 

- Hardware nonlinearities . . . control surface rate and deflection limits 

- Coded non linearities in digital control system 

• I h ® ob i® ctive to employ optimization to determine the excitation 
that produces the maximum gust loads on nonlinear aircraft 

- Matched filter theory for linear systems provides starting guess 


Figure 1 
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Schematic of Matched Filter Theory as Applied to Linear Systems 


Figure 2 contains a schematic of the steps necessary to generate maximum time-correlated gust loads for a 
linear system using Matched Filter Theory. The signal flow diagram is presented as two paths, the top path 
illustrates the generation of the system impulse response; the bottom path illustrates the generation of the 
maximum response of the system. 

In the top path, a prefilter (transfer function) representing the gust dynamics is excited by an impulse of unit 
strength to generate an intermediate gust impulse response which, in turn, is the excitation to the aircraft. 
Computationally, the time history of the response is carried out until the magnitude of the response dies out 
to a small fraction of the largest amplitude of the response. The response is normalized by its own root- 
mean-square (rms) value. This normalized response, reversed in time, is the matched excitation waveform 
for the output y. This becomes the input to the next part of the computational process. 

The bottom path illustrates how the maximum response of the system and the critical gust profile are 
obtained. The matched excitation waveform is applied to the same "known dynamics." The intermediate 
gust response is referred to as the critical gust profile. The final time histories are time-correlated gust loads. 
The maximum guaranteed response, ymax> * s equal to the rms of the impulse response and may be scaled to 
correspond to the gust intensity levels of interest. It should be mentioned that to obtain both the critical gust 
profile and the maximum response for a different output, a separate but similar analysis needs to be 
performed. 



von Karman Gust 
Impulse Response 


Time-Correlated Responses to 
Excitation Waveform Matched to y 



Figure 2 


223 



Optimization Scheme to Obtain Maximum Gust Loads for Nonlinear Systems 

The objective of the present research is to determine the excitation, with a prescribed rms, that results in a 
maximized peak gust load response while producing a gust profile of constant energy level (and thus a 
constant probability level) in an aircraft with a nonlinear element. This figure illustrates how an optimization 
algorithm may be employed to compute maximized gust loads and their corresponding critical gust profiles 
for nonlinear systems. The matched filter approach (as described in connection with figure 2) is used to 
provide an initial estimate of an excitation waveform for turbulence of a given intensity, shown in the upper- 
left comer of figure 3. The optimization scheme begins with the computation of the coefficients of a set of 
orthogonal functions in a series approximation to the waveform, normalized to a unit rms. The 
approximation to the excitation waveform is the input to the gust prefilter, whose output is an iterative gust 
profile. The gust profile then becomes the input to the nonlinear airplane model. The final output is a time 
history of the load quantity of interest. Note that the shaded area in the optimization loop is analogous to the 
"known dynamics" element in figure 2. 


The orthogonal approximating function coefficients, which are the design variables in the optimization 
scheme, are systematically varied by the optimizer until a maximum peak in the load response in obtained. 
The coefficients are constrained to produce a waveform approximation with a unit rms. Since the 
approximating functions are orthogonal, Parseval's Theorem allows the rms of the excitation waveform to be 
written simply as the sum of the squares of the coefficients. 
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Coefficient Generator Using Fourier Series Approximation 

Figure 4 presents an example of approximating the excitation waveform with the coefficient generator. 
Fourier series has been investigated as a candidate approximating function. The figure shows an initial 
waveform to be approximated, and two examples of Fourier series approximations. The second plot shows 
the resultant curve for 41 coefficients. The peak excitation is significantly underpredicted and there is a high 
frequency oscillation present during the latter portion of the time history. Using 401 coefficients to 
approximate the excitation sufficiently captures the curve's characteristics, as illustrated in the third plot. 


Original 

Waveform 



300 

200 

41 Term 100 
Approximation 



401 Term 



Figure 4 
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Block Diagram of Aircraft Control System with Nonlinear Element 

A simulation model of a drone aircraft was constructed to demonstrate an application of the present method. 
The model is based on a configuration used to design the gust alleviation control system as discussed in 
reference 3. Figure 5 shows the block diagram of the simulation model which includes the aeroelastic plant, 
the gust load alleviation control law, and the nonlinear control element. The shaded block to the left of the 
plant is the iterative gust profile input. The shaded block to the right of the plant is the wing root bending 
moment. The maximum absolute value of the wing root bending moment is the objective function. 

The plant itself is a linear, s-plane aeroelastic half-model consisting of 2 rigid body modes and 3 flexible 
modes. Unsteady aerodynamics were obtained using the doublet lattice method (reference 4). The plant 
model also includes the dynamics for the aileron and elevator control surface actuators. The two-input/two- 
output control law was obtained using a Linear Quadratic Gaussian design approach with the intent of 
reducing wing root bending moment. The controller uses the two control surfaces simultaneously. 

The original control system design did not contain any nonlinear elements. The nonlinear element defined in 
the figure is based on a spoiler-driven gust load alleviation system used on the Airbus A320 (reference 5). 
This nonlinearity is intended to simulate some of the important aspects of an actual system; these aspects 
include allowing motion only in one direction and preventing motion beyond a deflection limit. 

It should be added that wing bending moment response is dominated by the short period dynamics and is 
characterized by a large overshoot and a smaller undershoot. The objective of the load alleviation system is 
to reduce the overshoot load above a one g level and to ignore the undershoot loads below this level and the 
neutral load condition. The nonlinear element in the controller accomplishes this type compensated load 
reduction. 


Accelerometers 



Figure 5 


226 






Status 


Figure 6 outlines the status of the task. The entire scheme presented in figure 3 has been implemented. The 
nonlinear control system has been simulated with MATR1X X S YSTEM_BUILD (reference 6), which uses a 
high-level interpretive language and nonlinear functions that are built into the program. The simulation is run 
0I f a MicroVAX computer The nonlinear simulation with approximately 2000 time points takes about five 

minutes clock time to run. 

As indicated in figure 4, 401 terms in the Fourier series are necessaiy for an adequate representation of die 
initial excitation waveform. Since this is the number of terms used in the implementation, this task required 
the generation of 200 sine and 201 cosine waveforms for each objective function evaluation. 

The optimization is performed using a MATRIXx Optimization Module (reference 7) that also incorporates a 
high-level interpretive language. Gradients are generated from within this module using a finite difference 
method. This part of the computation is estimated to require 402 evaluations (one more than the number of 
design variables) of the objective functions which means performing the same number of simulations. 

Usine SYSTEM BUILD and the Optimization Modules, the optimizer was allowed to run for a day and a 
half clock time and had to be stopped. It was then decided to modify the method to allow more rapid solution 
of the problem Figure 7 outlines some of the problems encountered and the possible solutions. 


. Nonlinear control system implemented using MATRIX x 
SYSTEM_BUILD 

- Uses high-level interpretive language 

- Employs built-in nonlinear functions 

• Excitation waveform generator utilized Fourier series 

- 401 coefficients necessary for good initial waveform approximation 

- Waveform composed of 200 sines, 200 cosines, and 1 constant 

• MATRIXx Optimization Module used to maximize loads 

- Objective function is the peak wing root bending moment response 

- Uses equality constraint to maintain constant energy 

- Generates gradients using finite difference method 


Figure 6 
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Research Problems and Proposed Solutions 

Figure 7 presents the major research problems identified to date. The central issue is speed. To deal with 
this problem three areas are being investigated: simulation speed, number of design variables and optimizer 


The simulation constructed using SYSTEM_BUILD was run on a MicroVAX. Both the high-level 
interpretive language of SYSTEM_BUILD and the machine limitations of the MicroVAX contribute to the 
slowness of the simulation. To overcome this problem, a FORTRAN-based simulation needs to be 
generated using HYPER_BUILD (reference 8) and run on a faster computer using RemoteSim (reference 9). 

The coefficients used to generate the excitation waveform serve as the design variables in the optimization 
problem. Using Fourier senes to approximate the waveform requires an exceptionally large number of 
c 1 p e . lcl . en * s ;. Other orthogonal functions such as Chebychev polynomials are being investigated to determine 
their suitability for approximating the waveform. A reduction of the number of coefficients can be achieved 
by not approximating the discontinuous drop off to zero of the excitation waveform and explicitly setting 
waveform to zero at that point. r J B 


The speed of the optimization module is the third area for possible improvement. Increasing the speed could 
be accomplished by using a FORTRAN-based optimization module instead of the high-level interpretive 
language of MATRIXx. Since maintaining the equality constraint is a difficult task to achieve for most 
optimizers, the number of iterations through the optimizer loop could be reduced by reformulating the 
equality constraint as an inequality constraint. Gradients currently generated by finite differencing might also 
be generated analytically. The number of design variables could be reduced by using only the coefficients 
with the largest gradients. This would also produce a faster optimization. 


Problem: 

• SYSTEM_BUILD simulations too slow on VAX computers 
Solution: 

- Fortran-based simulation such as HYPER_BUILD needs to be generated 

- Simulation must be executed on faster machine 

Problem: 

• Exceptionally large number of Fourier coefficients are needed to generate the 
excitation waveform 

Solution: 

Use other orthogonal functions better suited for waveform approximation 

- Precomputed polynomial waveforms for later use 

Problem: 

• MATRIXx Optimization Module too slow 
Solution: 

- Incorporate a Fortran-based optimization module 

- Reformulate equality constraint as inequality constraint 

- Generate gradients analytically 

- Choose as design variables only those coefficients with the largest gradients 


Figure 7 
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PROCEDURE OUTLINE 


The purpose of this investigation is to develop an analytical method to study the 
vibration characteristics of piezoelectrically forced quartz plates. The procedure 
is schematically shown in Figure 1, and can be summarized as follows. The three 
dimensional governing equations of piezoelectricity, the constitutive equations and 
the strain-displacement relationships are used in deriving the final equations. For 
this purpose, a state vector consisting of stresses and displacements are chosen and 
the above equations are manipulated to obtain the projection of the derivative of 
the state vector with respect to the thickness coordinate on to the state vector 
itself. The solution to the state vector at any plane is then easily obtained in a 
closed form in terms of the state vector quantities at a reference plane. To 
simplify the analysis, simple thickness mode and plane strain approximations are 
used. 


THREE-DIMENSIONAL 
GOVERNING EQUATIONS 

( 


— 1 


EQUATIONS 


THICKNESS MODE 




AND PLANE STRAIN 

w 


STATE VECTOR 

APPROXIMATION 

W 

\ 

<4 

r 

SELECTION 


KINEMATIC 

RELATIONSHIP 


! DERIVATIVE OF STATE VECTOR 
WITH RESPECT TO THE 
THICKNESS AXIS 

f 

EXPLICIT EXPRESSION FOR 
STATE VECTOR AT ANY PL ANE 
IN TERMS OF THE STATE VECTOR 
AT A REFERENCE PLANE 


FU TURE WORK 


Figure 1 
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The governing equations of piezoelectricity consisting of the equations of motion 
and the charge equations of electrostatics are given by Equations (1) and (2). T e 
quantities , u^ and are the components of stress, mechanical and electrical 

displacements. The constitutive equations are presented in Equations (3) and (4), 
where C-hvi is the elastic stiffness, and efci, e^, and are respectively t e 

components of mechanical strain, piezoelectric strain constants , electric field and 
dielectric permitivity. The relationship between mechanical strain and displacement, 
and the relationship between electric field and electric potential are given in 
Equations (5) and (6) respectively. 


EQUATIONS OF MOTION 

a ij , j = p u i,tt 

CHARGE EQUATION OF ELECTROSTATICS 

Di.i = 0 

CONSTITUTIVE EQUATIONS 

^ij = c ijkl e kl " e kij E k 
D i = e ijk € jk + s ij E j 

Hj = °- 5 < u j,i + u i.j> 

Ei = 

Figure 2 


( 1 ) 


( 2 ) 


(3) 

(4) 

(5) 

( 6 ) 
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The plane xl-x3 is taken to be the plane of the plate, and the x2- direction is 
considered as the thickness coordinate. The simple thickness mode approximation, in 
which the various quantities are just functions of the thickness coordinate, is used 
m the analysis. Also, the system is considered to be under plane strain conditions. 
Invoking the above assumptions, and using the contracted notation given by Equation 
(7), the surviving system of equations are presented in Equations (8) through (10). 
Differentiating the last of the equations (10), using the third of Equation (8) and 
integrating the resulting equation twice, the expression for <f> is obtained 
(Equation 11), where A and B are constants of integration. A constant field does not 
produce any electric field, hence the constant B in Equation (11) is neglected. 
Substituting Equation (11) in Equation (10), the expressions for the non zero stress 
components are obtained, and are given in Equations 12 and 13. 


CONTRACTED NOTATION 


11, 22, 33, 23 OR 32, 31 OR 13, 12 OR 21 


I, 


(7) 


1. 2, 3,4, 5, 6 


INVOKING SIMPLE THICKNESS MODE AND PLANE 

a 6 , x2 " P U 1 , tt 
a 2 , x2 - P u 2 , tt 


STRAIN ASSUMPTIONS 


( 8 ) 


°2,x2 - 0 


€ 2 “ u 2>x2 
e 6 - 1/2 u 1|X2 

E 2 " -<*.x2 


(9) 


a 2 ~ V2 C 26 u l x2 + C 22 u 2 x2 + e 22 ^ x2 

a S “ c 66 u l,x2 + c 62 u 2,x2 + e 26 ^,x2 

°2 - V2 e 26 u l x2 + e 22 u 2 x2 - S 22 * >x2 


( 10 ) 


<t> - ( 1/2 e 2 g +e 22 u 2 )/S 22 + Ax2 + B (11) 

a 2 ” a 26 U 1 , x2 + a 22 u 2 , x2 + e 22 A (1 2 ) 

a 6 ~ a 66 U 1 , x2 + a 62 u 2 , x2 + e 26 A (13) 


Figure 3 


234 



A state vector (V) defined by Equation (14) is chosen. The derivatives of the state 
vector with respect to x2 is obtained from Equations (8) through (13) and the 
resulting expressions are given in Equations (17) through (20). The elements of the 
matrices Bl, B2 and B3 are made up of the material constants and derivatives with 
respect to x2 and time (t). 


STATE VECTOR 


(V) - [ ( V L ) T 

[VI) ” [ U]_ &2 

(V 2 ) T ] t 

] T ; { V 2 ) = [ a 6 u 2 ] T 

(14) 

(Vl} )X 2 - I Bi 

] { V 2 ) - A {b x ) 

(15) 

{V 2 } iX 2 - [ b 2 

] (V!) - A {b 2 } 

(16) 

m, x2 = [ B 3 

] (V) - A (b 3 ) 

(17) 



l/ a 66 

-( a 62/ a 66) d/dx 2 


e 26/ a 66 

Bl = 

0 

p d 2 /dt 2 

bl - 

0 


p d 2 /dt 2 0 


0 


b2 - 


-( a 26/ a 22) 3/3 x 2 1/ a 22 


e 22/ a 22 


(18) 


(19) 


2 

a 62 “ c 62 +e 26 e 22/ s 22 < a 66 “ V 2 < c 66 + e 26 / s 22> > 

2 

a 26 = V 2 { C 26 + e 22 e 26/ s 22 ) I a 22 “ c 22 + e 22 / s 22 


( 20 ) 


Figure 4 
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A solution to the differential equation given in Equation (11) can be easily 
obtained and is given in Equation (21), where {VO} is the state vector evaluated at 
x2-0 . The analyst has the flexibility of choosing any plane as the appropriate 
reference plane. The exponential term in Equation (21) can be expressed in an 
infinite series , and the powers of the matrix B3 can conveniently grouped as shown 
in Equation (23). 


v - e B3 x2 {VO} + B3 ’ 1 A b3 

(VO) = {V} x2=0 


e B3x2 


(VO) 


I + B3 x2 + (B3 x2) 2 /2!+(B3 x2) 3 /3!+ 


(VO} (22) 


p 

0 1 

b3 3 . r 0 

B1Q -1 

- 0 

Q J 

L B2P 

o J 


^ P On 

5 r 0 

B1Q 2 - 

A 

B3 3 — 


0 Q 2 J 

L B2P 2 

0 - 


(23) 


[ P ] = [Bl] [B2] 
IQ]- [B2 ] [Bl] 


Figure 5 
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Using the expressions given in Equation (23) , the infinite series expansion for the 
exponential term can be conveniently grouped as shown in Equation (24). The elements 
present in Equation (24) can be recognized as a convergent series. The resultant 
expression is given by Equation (25). Substituting this expression in Equation (21), 
the final equation for the state vector at any reference plane in terms of the state 
vector at a reference plane is obtained (Equation 26) . 



0 B1Q 2 
- B2P 2 0 


x2 5 /5 ! + 


{VO) 


(24) 


cosh (x27P) B1 70 sinh (x27Q) 

B27P sinh (x2./P) cosh (x27Q) 


(VO) 


(25) 


- [R] {VO) 

(V) = [R] (VO) + A [ B3 ] " 1 {b3 } 


(26) 


Figure 6 
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Symbols and Abbreviations 


A, B Constants of integration 
a ij Constant coefficients 

c ijkl Elastic stiffness 

&i Components of electric displacement 

e kij Components of piezoelectric strain constant 

E^. Components of electric field 

Sij Components of dielectric permittivity 

u i Components of mechanical displacement 

Components of strain 
<t> Electric potential 

p Mass density 

j Stress components 
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Introduction 


For more than two decades, viscoelastic materials have been commonly used as a 
passive damping source in a variety of structures because of their high 
material loss factors. In most of the applications, viscoelastic materials 
are used either in series with or parallel to the structural load path. The 
latter is also known as the constrained- layer damping treatment^ » ^ . The 
advantage of the constrained- layer damping treatment is that it can be 
incorporated without loss in structural integrity, namely, stiffness and 
strength. However, the disadvantages are that (1) it is not the most 
effective use of the viscoelastic material when compared with the series-type 
application, and (2) weight penalty from the stiff constraining layer 
requirement can be excessive. To overcome the disadvantages of the 
constrained- layer damping treatment, a new approach for using viscoelastic 
material in an axial-type structural components, e.g., truss members, was 
studied in this investigation. 
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Elastic Tailoring in Composite Structures 


It is well known that, with the properly arranged orientation sequence in 
layup, composite structure can exhibit various types of deformation coupling 
when subjected to loading. In certain applications, such anisotropic behavior 
can be tailored to benefit specific needs. For example, the bending/ twisting 
coupling has been extensively studied for the purposes of aeroelastic 
tailoring-^ » ^ . The application of extension/twisting deformation coupling to 
the constrained- layer damping treatment was explored in Ref. 2. In Refs. 5 
and 6, a new approach of applying extension/twisting deformation coupling to 
damping treatment was proposed for the axial -type truss member. In this 
approach, the viscoelastic material is embedded in a structural member made of 
fiber reinforced composite material. By a judicious tailoring, the structural 
member can exhibit the extension/twisting deformation coupling such that the 
viscoelastic material is sheared in twisting while the structural member 
undergoes an axial deformation. 
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New Material Concept: 

However, the difficulty with this new approach is that it requires a built-in 
twisting freedom in the truss member. In reality, such added design 
complexity is usually forbidden. To avoid such undesirable requirement, a new 
material concept of using saw-toothed (or V-shaped) fiber in reinforced 
composite was conceived in Ref. 7. With the V-shaped fiber, a truss member is 
allowed to undergo twisting deformation at knee-points while its both ends 
remain fixed. Damping performance was studied on a plane strain model as 
shown below. The resulting shear strain distribution in the viscoelastic 
material is a hyperbolic sine function along the member axis which is similar 
to the result of constrained- layer damping treatment 2 . The member loss factor 
is estimated from the expression of complex stiffness as 

_ tanh Imaginary (K m ) 

K m - K - ( 1 + fi- ) iy m = 

■IP Real (K m ) 

where 


Al6 2 

V = 

a 11 a 66 - a 16 2 


P = L 2 • 


2G 2 


c 2 


A 11 

a 11 a 66- a 16 2 


Examining the above expression indicates that the member loss factor is only a 
unction of two parameters. One is the extension/twisting coupling 
coefficient, /i, which is a function of the composite material properties and its 
iayiip. The other is a combined geometry/material parameter, p. 
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Analysis Model for Waved Fiber Reinforced Composites 


The idealized V-shaped fiber is useful in performance trade study. In 
practice, however, the V-shaped fiber is not really feasible because of its 
sharp corners. In this investigation, the sine-waved fiber reinforced 
composite is analyzed for this new damping treatment approach. The analysis 
model of a truss member made of iso-phased sine-waved fiber reinforced 
composite is shown below. In this design, the fiber orientation is 
antisymmetric with respect to the viscoelastic material (VEM) layer, i.e., 
[+0 n /VEM/-0 n ] layup, such that the twisting deformation in the viscoelastic 
material is maximized under axial load. Because of the continuously varying 
fiber orientation, the truss member's elastic properties are varying along the 
member axis. In this study, a concept of using equivalent homogeneous model 
with effective elastic properties is proposed to evaluate the member's damping 
performance. The effective anisotropic properties are estimated in average 
sense over one quarter of the fiber's spatial wavelength. 




EFFECTIVE 
ANISOTROPIC 
PROPERTIES: 
Sjj, Qjj, Ay 
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Spatial dependent strain-stress relation: 


(«> = [S(0)J (a) 

where the S-y's are the elemental elastic compliance and the tangent fiber 
orientation, 6, is given by 

2na 

0 = tan • cos(27rx/A) ] 

A 

The effective compliance is 
1 rA/4 

s ij = * Sjj (x) * dx 








Given the material loss factor, rj vem -1.0, functional dependency of the 
member loss factor on the extension/twisting coupling coefficient, fi, and 
parameter y3 are illustrated in the following example. It is interesting to 
note that the /3 0 pt is not ver y sensitive to the variation in n. For the 
example of HMS/3501-6, /x G pt = 1.2 at r = 18°, the maximum loss factor 
attainable is about 17%. 



3 

2 



Summary 


o New material concepts, i.e., V-shaped and sine-waved fiber 
composite materials, were investigated for the new damping 
axial-type structural members. 


reinforced 
treatment in 


The underlying mechanism of an elastically tailored damping member depends 
on he extension/ twisting deformation coupling in composite materials As a 
result, the embedded viscoelastic material is sheared in twisting when the 
member undergoes axial motion. 6 


o 


Shear strain distribution in the viscoelastic material is 
the extension/twisting coupled damping treatment and the 
damping treatment. 


similar between 
constrained- layer 


o A concept of using an equivalent homogeneous model with effective 

anisotropic properties was proposed to evaluate damping performance of 
members made of iso-phased sine-waved fiber reinforced composite material 


o Numerical examples show that the sine waved- fiber reduces the degree of 
extension/ twist mg deformation coupling as compared with the V-shaped fibe 
reinforced composite material. However, its effect on the 0 parameter is 
less critical because the p parameter can be optimized through other 
geometric parameters. 


r 


o With the optimally selected geometric and material parameters the 

?rom^n n S w t 0 r K° f ^ elasticall y tailored damping member ranges 

om 10-25% which is about in the same performance range of the damping 

member with constrained- layer damping treatment. The major advantage of the 
elastically tailored damping member is that there is no additional weight 
penalty such as the constraining layer of the constrained- layer damping 
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OBJECTIVES 


This work-in-progress presentation describes an ongoing research activity at the NASA Langley 
Research Center to develop analytical methods for the prediction of aerothermoelastic stability of 
hypersonic aircraft including active control systems. The objectives of this research include application of 
aerothermal loads to the structural finite element model, determination of the thermal effects on flutter, and 
assessment of active controls technology applied to overcome any potential adverse aeroelastic stability or 
response problems due to aerodynamic heating- namely flutter suppression and ride quality improvement 
K>r this study a generic hypersonic aircraft configuration was selected which incorporates wing flaps 
ailerons and all-moveable fins to be used for active control purposes. The active control systems would 
use onboard sensors in a feedback loop through the aircraft flight control computers to move the surfaces 
or improved structural dynamic response as the aircraft encounters atmospheric turbulence. 


• Construct a Generic Hypersonic Vehicle to Use in 
Studies 


Performing Analytical 


• Develop and Analyze Aeroelastic Models Incorporating the Effects of 
Aerodynamic Heating 


• Apply Active Controls to Compensate 
for Degraded Dynamic Responses 


• Flutter Suppression System 

• Ride Qualities Augmentation 
System 
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HYPERSONIC ANALYSIS AND DESIGN APPROACH 

The current aeroservothermoelastic (ASTE) analysis and design capability is outlined schematically 
below. The method consists of three primary steps; 1) the determination of thermal loads acting on the 
structure due to aerodynamic heating, 2) the development of hot and cold aeroelastic mathematical models 
for flutter analysis including the computation of unsteady aerodynamic forces acting on the structure, and 
3) the design, analysis, and simulation of active control laws. 
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APAS AEROTHERMODYNAMIC MODEL 

c T|fHy^ onic Arbitrary Body Program (HABP) of the Aerodynamic Preliminary Analysis 
bystem 1 1 J (APAS) was used to model the generic hypersonic aircraft configuration and obtain steady- 
state aerodynamic forces and heat loads. For a given flight condition (angle-of-attack and control surface 
aettection), the HABP module was used to compute aerodynamic lift and moment coefficients and 
aerodynamic center location, as well as the radiation equilibrium wall temperatures on the vehicle. The 
aerodynarmc results were used to calibrate the unsteady aerodynamic force calculations by comparison of 
pitching moment coefficient and aerodynamic center location. The unsteady aerodynamic force models 
were then modified to yield compatible results. The radiation equilibrium wall temperatures were used 
directly as heat loads m the finite element structural model to determine structural stiffness changes caused 
by thermal stresses and material property changes. 



Aerodynamic Preliminary Analysis System (APAS) Hypersonic 
Arbitrary Body Program (HABP) module used for steady-state 
aerodynamic calculations 

• Radiation equilibrium wall temperatures 

• Lift and moment coefficients, aerodynamic center locations 

Results used to 

• Provide heat loads for thermal structural analysis 

• Calibrate unsteady aerodynamic codes 
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FINITE ELEMENT MODEL 


A conventional structural concept was used for the generic aircraft configuration of this study [2]. 
The fuselage was modeled as an elliptical cross section (width/height ratio 2) consisting of stiffened ring 
and skin construction. The low-aspect wings were modeled as fully attached to the fuselage consisting of 
spars, ribs, and skins. The wing leading edge sweep is 70 deg. and the wing section is a 3% circular arc 
airfoil. A body weight fraction, defined as the weight of the structural material contributing to stiffness 
divided by gross takeoff weight, of 8.6% was used to determine the required structural mass. Material 
properties consistent with titanium aluminide were assumed for all structural elements. The wing flaps, 
ailerons, and all movable fin were modeled separately and attached to the fuselage/wing model by spring 
stiffness elements modeling actuator stiffness characteristics. The Engineering Analysis Language [3] 
(EAL) structural analysis code was used to compute hot and cold vibration mode frequencies and mode 
shapes. The visual appearance and overall character of the mode shapes did not change with variations in 
temperature, although significant changes did occur in frequencies. Heating effects decreased the 
frequencies by thirteen to twenty percent. 



NATURAL FREQUENCIES (Hz) 


MODE COLD HOT 


3.01 

2.43 

4.02 

3.48 

7.06 

5.67 

7.70 

6.56 

9.47 

7.63 

10.96 

8.84 
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UNSTEADY AERODYNAMICS - LESSONS LEARNED 

Significant problems were encountered in computing valid unsteady aerodynamic forces for use in 
aeroelastic stability analyses in both subsonic and supersonic flight regimes. For the subsonic case, two 
versions of the Doublet Lattice Method [4] (DLM) aerodynamic panel code were used, as was a Kemal 
Function Method [5] (KFM) code. In the case of the DLM, the two versions were inconsistent in force 
results (both magnitude and phase). This was attributed to nonconvergence of the DLM due to insufficient 
numbers of aerodynamic boxes. The minimum number of required boxes was later estimated to be on the 
order of 675, far exceeding reasonable computational cost. Subsonic flutter boundary predictions using 
the KFM code were erratic, showing wide oscillations in flutter dynamic pressure for small subsonic 
variations in Mach number. For the supersonic case, the MSC/NASTRAN [6] Mach Box [7] and Piston 
Theory [8] methods were tried. It was found that the Mach Box result would not compare with analytical 
solutions for simple check cases. The Piston Theory method was found to be restricted to rigid chords, 
typically valid for high aspect ratio wings which are very stiff chordwise, and did not include airfoil 
thickness effects. Two new second-order Piston Theory codes including thickness, camber and 
chordwise bending effects were written, one in EAL and one in FORTRAN, both taking advantage of an 
ousting aero/structure interface [9]. The FORTRAN version aerodynamic force results were ultimately 
used for flutter analyses because of consistency with the earlier APAS steady-state results. 


SUBSONIC: 

Doublet Lattice 

• Inconsistent between code versions (ISAC and NASTRAN) 

• Estimated 675 aerodynamic boxes required for convergence 

• Exceeds inhouse code capability, very expensive in NASTRAN 

Kernel Function 

• Erratic flutter boundary predictions 

SUPERSONIC: 

Mach Box 

• NASTRAN results do not agree with analytical solutions for 
simple cases 

Piston Theory 

• NASTRAN mode! limited to (a few) rigid chord panels 

SUPERSONIC SOLUTION: Write a new piston theory code 

• Linked to ISAC aero/structure interface to model nonrigid chords and 
thickness 
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PISTON THEORY AERODYNAMIC IMPLEMENTATION 

At sufficiently high Mach numbers "local" wave theory is a good approximation to the unsteady 
aerodynamics. The local pressure is related to the normal free stream velocity in a similar manner as the 
pressure in a one-dimensional piston chamber is related to the velocity of the piston. A local, linearized 
pressure equation is represented by the equation shown in the figure. The various aircraft surfaces were 
represented by trapezoidal panels similar to the one indicated in the figure. The normal velocities over the 
surfaces were computed using surface spline interpolation with the normal velocities located at the center 
of each trapezoidal panel. The point forces subsequently created by the piston theory pressures were also 
concentrated at the center of each panel. The generalized aerodynamic force for each mode was generated 
by summing these point forces, weighted by the interpolated mode shapes, over the aircraft surfaces. The 

additional symbols used in the figure are defined as follows: Ap, pressure difference between upper and 

lower surfaces; p, density; a, the local speed of sound; y, ratio of specific heats; Z, the relation describing 
the contour, or the thickness of the vehicle component; z, the displacement of the discrete point; and V, the 
freestream velocity. 



(U) Linearized, second-order equations including thickness effects 

Ap(x,y,t) = -2pa 1 + G ;! Z(x ' y> {' V 3x + 

r M 4 (v+1)-4P 2 
2(3 

(U) z(x,y,t) calculated at discrete points using surface spline interpolation 
of mode shape data 

(U) Generalized aerodynamic force for each mode computed by numerical 
integration of the pressure over the surface 


; p = VM 2 
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AERODYNAMIC MODELING INFLUENCE ON FLUTTER DYNAMIC PRESSURE 

The relative importance of the aerodynamic influence of various vehicle components to flutter was 
evaluated. The importance of the inclusion of a modeling effect was measured by the percent change in 
flutter dynamic pressure. Four of the effects examined were found to be significant. The baseline 
analysis contained a restrained flat plate representation of the wing and used the first six flexible modes. 
Introducing the rigid body plunge and pitch modes into the analysis increased the flutter dynamic pressure 
by ten percent. Because the structural frequencies are very low, the structural modes are influenced by the 
short period mode. For this configuration, the rigid body motion helps to dissipate the system's energy 
into the airstream, thus inhibiting flutter. The addition of a flat plate representation of the fuselage 
decreases the flutter dynamic pressure by ten percent. Fuselage motion dominates the first and third 
flexible modes; including this motion in the analysis increases the aerodynamic force input and encourages 
flutter. Addition of an aerodynamic representation of the vertical fin further decreases the flutter q by ten 
percent for the same basic reason. The fifth flexible mode has significant vertical fin contributions, 
milking it important to the analysis. The remaining changes to the aerodynamic model are inclusion of the 
thickness effects of the wing and the fuselage, both of which cause an increase in the flutter dynamic 
pressure. The wing contour effects changed the flutter value by ten percent, while the fuselage contour 
effects changed it by only two percent. The final model used for analysis and design incorporated all of 
the above effects except for fuselage thickness. 


Basic Model; 

Flat Plate Clipped Delta Wing with 70 Degree Leading Edqe Sweep 
First 6 Flexible Modes 


CHANGE TO MODEL 

■ Inclusion of Rigid Body Pitch & Plunge Modes 

■ Addition of Flat Plate Fuselage 

■ Addition of Wing Thickness Effects 

■ Inclusion of Flat Plate Vertical Fin 

■ Addition of Fuselage Thickness Effects 
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ALTITUDE ROOT LOCUS 


This figure shows the variation of the eigenvalues associated with the rigid body and first four 
flexible modes, as altitude is varied. The model used for this typical root locus was the heated structure at 
Mach 2.0. The arrows indicate decreasing altitude. This analysis was performed at a matched point, so 
both the density and velocity changed with altitude. Flutter is determined as the cross-over of the 
imaginary axis. From the figure it is seen that the eigenvalue associated with the third flexible mode 
moves into the right half plane at a frequency of 43.5 radians per second. 


MACH 2 HOT 


IMAG. PART 
(1 / SEC) 



REAL PART 
(1 / SEC) 
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EFFECTS OF MACH NUMBER AND HEAT ON SHORT PERIOD MODE DYNAMICS 


The short period mode dynamics are influenced by both the structural properties and by the 
aerodynamics. The figure shows the changes in short period behavior incurred due to the effects of 
structural flexibility, heating and Mach number. The six curves represent the trace of the short period 
eigenvalue in the complex plane as the altitude is varied ( altitude decreases from right to left along each 
curve). 

In an aeroelastic system, the roots of any one mode are influenced by the other modes near it. 
Because the structural frequencies for these configurations are low, in the neighborhood of the rigid body 
frequency, it is anticipated that they would exhibit a large degree of influence over the short period mode. 
This influence can be seen by examining the roots of the rigid vehicle versus the eigenvalues after the 
effects of flexibility have been included. The figure indicates that including flexibility tends to have a 
destabilizing effect . The effects due to the aerodynamic heating can be seen by comparing the hot and 
cold data for the same Mach number. At either Mach number, the destabilizing effect of the heating is seen 
as the roots for the hot data fall further to the right in the s-plane than those corresponding to the cold data. 

To determine the effects of the aerodynamics, the curves for the rigid data, the hot data and the cold 
data must be examined separately. It is seen that as the Mach number is increased, the short period 
frequency is increased and the damping is decreased. Thus, increasing Mach number also has a 
destabilizing effect on the short period dynamics. 

Comparing the curves in these three ways shows clearly that the Mach number has a much larger 
influence than either the flexibility or the heating. 


Short period root locations at various altitudes 
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-0 COLD FLEXIBLE 

HOT FLEXIBLE 

MACH 2 

RIGID 

COLD FLEXIBLE 
HOT FLEXIBLE 
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Real Part (1/sec) 


Imag. Part 
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HOT / COLD FLUTTER RESULTS 


The flutter characteristics are presented as a set of curves showing regions of instability. e 
flutter boundaries illustrate the destabilizing effects of both heating and Mach number. The region ow 
either curve represents the region for flutter-free flight. As the Mach number is increased, die flutte 
dynamic pressure is increased. Heating lowers the flutter boundary over the entire range o ac 
numbers, indicating that there will be an instability at lower dynamic pressures. 


FLUTTER BOUNDARIES FOR 
HEATED AND UNHEATED 
VEHICLES 


1.50 
1.25 

1.00 

q / q ref o.75 
0.50 
0.25 

0.00 

12 3 4 

Mach Number 



259 



CONCLUSIONS / FUTURE PLANS 


, rin(T aer0therm0 ! lastic method has been developed. The thermal loads due to aerodynamic 

r the r Tu e n lemen i ^^y 515 - A PP licad on of the aerothermal effects reduced 
structural frequencies and lowered the flutter boundary. The flutter was found to be influenced by all 

wing flmter ° f ^ VehlC e; anal y ses of hypersonic configurations must consider aircraft flutter versus 

,H^ UtUre w ° rk in this aica wi ? concentrate on control law design and closed loop analysis. Plans 

to thal of^e S coH S p^ reSS1( ; n syStem whlch wil1 raise flutter boundary of the heated vehicle up 

S It!!! , the cold vehicle. Ride qualities improvement will also be a focus in the control law design phase 

iminveH Je S' ^ dd i!i 10n w^ in u ear unstead y aerodynamic codes will continue to be evaluated and P 

coE^s .?hyp e ™„i?SSis r Kd and U " hea,ed VeWCleS Wi " be deflned fom ,he subso ™ ni 8 h « 


CONCLUSIONS 

• Aerothermoelastic Analysis Method Developed 

Aerothermal Loads Incorporated into Finite Element Analysis 

• Reduced Structural Frequencies 

• Lowered Flutter Boundary 

• Flutter influenced by all vehicle components 

• Must Consider Aircraft Flutter Instead of Wing Flutter 


PLANS 


Further Evaluate and improve Linear Unsteady Aerodynamics Codes 


• Define Flutter Boundaries for Hot and Cold 
Hypersonic Speeds 


Vehicles from 


Subsonic to 


Hyttnte A?rcraft r °' S *° De,i " e Technol °gy Be "efits for 
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rro'foi BUL,L STATISTICAL ANALYSIS OF FIELD INSPECTION AND AIRCRAFT 
USAGE DATA HAS BEEN USED TO PREDICT THE RISK OF STRUCTURAL FAILURE 


We have described in previous work (ref. 1 and 2) the use of damage 
tolerance analysis and Weibull statistical analysis in the 
assessment of structural risk. The interference of the failure 
distribution and the aircraft life distribution is computed to 
determine the risk of structural failure. Information from any 
number of aircraft from different bases can be combined to give a 
projection of the risk associated with continued operation at the 
same or modified usage levels. 


Three parameter Weibull distributions are determined from the flight 
usage data and the failure information obtained from field 
inspection of the aircraft. In the present analysis, deterministic 
flaw growth analysis is used to project the failure distributions 
from inspection data. 



Figure 1 
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DETERMINATION OF FAILURE DISTRIBUTION FROM FIELD 
SERVICE INSPECTION DATA 

Inspection data is reported for each critical point in the aircraft. 
The data will indicate either a crack of a specific size or no 
crack. The crack length may be either less than, equal to, or 
greater than critical size for that location. 

Non-critical length cracks are projected to failure using the crack 
growth characteristics for that location to find the life when it 
will be at critical length. Greater-than-crit ical length cracks are 
projected back to determine the life at failure, that is, when it 
was at critical length. The same process is used as in the case of a 
non-critical crack except that the projection goes the other 
direction. These points, along with the critical length cracks are 
used to determine the failure distribution. 

To be able to use data from different aircraft to build a common 
failure distribution, a consistent life variable must be used. 
Aircraft life varies with the severity of the usage, therefore the 
number of flight hours for a particular aircraft must be modified by 
its usage factor to obtain a normalized life which can be compared 
with that from other aircraft. 


CRACK 

SIZE 
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USAGE FACTOR ALLOWS THE COMPARISON OF DATA FROM DIFFERENT AIRCRAFT 

The aircraft is designed to a baseline or design spectrum. This is 
determined from the design mission requirements for the aircraft 
The actual usage of the aircraft will vary greatly depending upon 
where the aircraft is based when it enters service. Some bases fly 
many more benign flights and others fly more severe flights than the 
baseline. For flight hours to be compared from one aircraft to 
another, they must be related to the same severity level or no 
direct comparison is possible. The usage factor is used to adjust 
the actual number of flight hours for the difference between the 
baseline usage and the actual usage of the aircraft. This method has 
been shown (ref. 3) to accurately account for the effect that usage 
has on the crack growth characteristics. The usage factor is the 

th^baseline ° f ^ SirCraft f ° r the »»•«• to 


L. L 

UF S = l^(> 10) UF m = ^-(< 1.0) 

s L m 
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FIELD DATA IS USED TO DETERMINE THE THREE-PARAMETER 

WEIBULL DISTRIBUTION 

Data from field inspect ions are used to determine the failure and 
life characteristics of the aircraft under consideration. The 
distribution of current lives is found from the number of hours 
(adjusted by usage) recorded for each aircraft. The failure 
distribution is found from the set of lives associated with the 
critical crack lengths. Again, the lives must be adjusted for the 
difference in usage. 

Linear regression is used to determine the best 3-parameter Weibull 
fit to the data. The median ranks are determined for the failed 
points and take into account the effects of the suspended items 
(non-cracked aircraft) on the rank values. The minimum expected life 
is found from a search process which determines what minimum life 
value gives the best straight line fit to the data. 

The difficulty with this process is twofold. First, there are 
generally only a few cracked parts from which you want to construct 
the failure distribution. The accuracy of the distribution so 
computed can be questioned. Second, the growing, or projecting, 
process assumes that the crack growth characteristics are 
deterministic . 


Percent 

Failured 



Figure 4 
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MULTIPLE FAILURE MODES ARE SOMETIMES PRESENT 

failures will sometimes result from several phenomenon. 

Manufacturing or material defects can precipitate early failures. 
These will generally occur well before the normal service failures. 
These failures are of interest, but it is important to separate this 
behavior from the normal service behavior for fleet management 
purposes. In addition, it is improper to attempt to fit a Weibull 

distribution to the combined data set since it does not correctly 

character ize either behavior pattern. The data set must be pruned to 

include only the long-term effects of the normal service life if an 

accurate picture of the failure rate and risk are desired. Generally 
the bulk of the data will be in this set, with the early failures 
being few in number. 

Similarily, if one wants to concentrate on short-term failures, the 
data must be pruned of other failure modes. Plotting all data, as 
shown in this chart, can help identify when more than one failure is 
represented in the data. 
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INITIAL INSPECTION DATA FOR 158 AIRCRAFT SHOWS 6 FAILURES 

Inspection of 158 fighter aircraft revealed the existence of 6 
aircraft with cracks of critical length at a point of concern on the 
vertical tail. Computation of the Weibull distribution shows that 
the data fits the curve fairly well, exhibiting a 0.97 correlation 
coefficient . 

Closer examination of the data points indicates that perhaps there 
are two failure modes present. The first failure at 770 hours seems 
to be isolated from the remaining five points. 



Figure 6 




CUMULATIVE PROBABILITY OF FAILURE FOR ORIGINAL DATA 

The cumulative probability of failure for the original data set 
ontaining six failures is shown. Included on the plot is the 90% 

5ecfsfnn Ce J and - The Confidence band is very important to the 
decision making process since frequently (as in this case) there are 

decision**** failures from whlch the fleet commander must reach a 

The confidence bands were computed using two different methods The 
f!ve and ninety five percent ranks were computed and fi? with i 

d Jf tributl0 £ alon 6 with the median ranks. This method 
provides the range for all three Weibull parameters; however the 
computation of the ranks and the curve-f itting procedure result in a 
substantiai computation time. The second method utilized the^ 
distribution to compute the confidence band for the linear 
regression parameters for the curvefit to the median ranks This 

Weibull locatio faster I how ever, we obtain no information for the 
loiatiL 1 ^ ? Parameter. This is a significant loss because the 
*5? r re P resents the failure free operating period The 

is fen to ou?ielgh g tSn a nss° n£1,JenCe li " it ‘ f ° r th * “''ail.ble data 


5% Curve 

Median Curve 

95% Curve 
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A SUBSEQUENT INSPECTION INCREASED THE DATA SET TO 181 AIRCRAFT 

WITH 12 FAILURES 

Subsequent inspection data increased the sample to 181 aircraft 
containing 12 aircraft with failures. Again this information was 
plotted and Weibull distributions determined for the median, five 
percent, and ninety five percent rank points. These curves are shown 
along with the result obtained by computing the confidence bands for 
the linear regression parameters. The two methods compare well, 
except at the lower end where the variation in the location 
parameter is felt more strongly. 
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CUMULATIVE PROBABILITY OF FAILURE FOR SECOND DATA SET 

The cumulative probability of failure for the second data set 
containing twelve failures is shown. Included on the plot is the 90% 
confidence band. 

The 90% confidence band is much smaller than that with only six data 
points, especially at the high probability of failure, indicating 
that the data set now represents the actual behavior of the failure 
mechanism to a much higher degree than the original data set. The 
influence of the early failure has been reduced by the new data 
points, many of which fell between the first failure at 770 hours 
and the second failure at 1035 hours in the original set of data. 


Median Curve 

■ 5% Curve 

— 95% Curve 



Flight Hours 


Figure 9 
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AREAS OF CONTINUING EFFORT 


We are continuing our effort in several areas. We will implement a 
Maximum Likelihood Estimation (MLE) process to determine the Weibull 
parameters. An iterative procedure is required; however, our 
experience with the MLE process for two-parameter Weibull curvefits 
indicates that convergence is very rapid. The linear regression 
process we are currently using weighs all the points equally in 
their effect on the regression line, whereas the MLE process weighs 
the analysis toward the bulk of the data. 

The process of projecting cracks to their critical level is 
accomplished deterministically from the crack growth curve. The 
crack growth process is, in fact, a random process and thus there is 
some uncertainty associated with the actual lives at failure. 
Inspection data is also treated deterministically. Nondestructive 
Evaluation (NDE) techniques have some uncertainty associated with 
their ability to detect flaws. The uncertainty, or randomness, of 
these two phenomena should be included. This uncertainty is best 
addressed using a Monte Carlo technique at the cost of some 
additional computation time. The advantage is that we will receive a 
better picture of the actual risk. 

Our current process does not account for the repair of cracked parts 
and the return of the aircraft to service. We are looking to Renewal 
Analysis techniques to provide an assessment of such repairs. 
Repaired aircraft are of particular interest to fleet commanders in 
planning allocation of resources and logistic needs and to project 
the maintenance and repair actions required with continued fleet 
usage . 


Maximum Likelihood Estimation of Weibull Parameters 

Monte Carlo Simulation Will Allow: 

- Random Crack Growth Characteristics 

- Uncertainty in Inspection Data 

Renewal Analysis to Account for Replacement of Failed 
Components 

Figure 10 
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PROBLEM DEFINITION AND CORRESPONDING RESEARCH 

The research is focused on automating the evaluation of comnlev 
structural systems , whether for the design of a neS system or the 

teohn? 1S v. an ® Xisting one ' b y developing new structural analysis 

ba J ed ° n qualitative reasoning. The problem is to identify 

and 2? the 11 the Requirements for the automation of design 

escribe the design and cognition components of the research. 

Design has received special attention in cognitive science ber*n<= e il- 
ls now identified as a problem solving activity that iJ 

various reasoning processes involving different kinds of knnwipHno • 
ways which vary from one context to anothe? ?L oSjec^ve Ts ^unUy 
all the various types of knowledge under one framework of cogni?ion * 

This presentation focuses on the cognitive science framework that tto 
knowledge and engineering reasoning in design. ^ 


RESEARCH: 


EVALUATE THE AUTOMATION OF COMPLEX STRUCTURES 
APPLICATION TO DESIGN AND ANALYSIS 


PROBLEM TO SOLVE: 


OBJECTIVE: 


UNDERSTANDING QUALITATIVE REASONING 


INTEGRATED TOOL FOR DESIGN AND RISK 
ASSESSMENT 


PRESENTATION: 


COGNITIVE ASPECT OF DESIGN 


DESCRIPTION: 


APPLICATION: 


THE MULTIPLE LAYER SEMANTIC NET — 
DESCRIPTION OF THE TYPE OF KNOWLEDGE IT 
SHOULD HANDLE 


ENGINEERING KNOWLEDGE 


DESIGN KNOWLEDGE 

Figure 1 
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KNOWLEDGE: PROCESSING OF CONCEPTS 


The common denominator among diverse entities such as an overall 
complex system, a component or a sub-assembly of that system, and the 
design and evaluation processes themselves, is that they can all be 
represented by formal concepts which, being associated with the human 
mind, can fundamentally encapsulate models of the reality that 
surrounds us [2] (percepts and icons). Concepts are organized in 
conceptual graphs, semantic nets, and schema or prototypes. Procedures 
can also be represented in semantic nets [7]. 

Different design reasoning procedures could be represented in various 
refinements of the same higher-order semantic net which corresponds, 
at the highest qualitative level, to deriving the structure for a 
device such that the device can meet a specific function. 


SYSTEM 1 


SCHEMAS / PROTOTYPES 


SEMANTIC NETS 


CONCEPTUAL GRAPHS 


CONCEPTS 

r 7 

PERCEPTS AND ICONS 

L J 


SYSTEM 4 


SYSTEM 2 


SYSTEM 3 
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CONCEPTS AS SEMIOTIC PARADIGMS 

Each concept associated with a cognitive process has three fundamental 
components: A semantic component to describe its function (what it is 
£ or >' a syntactic component to describe its structure (how it is out 

con?SJ r } # h a J d 3 component to describe how it relates to its 

ontext (what are its behavior and the context in which it is used) 
Pearson [3] attributes these components to cognitive systems and calls 
such concepts semiotic paradigms. 

Physical symbol system [4] and the connection models [5] have the 
f J heir Paradigms, but vary by the emphasis on the 
level of representational abstraction at which they are described. 

models of a device and the corresponding knowledge can be 
- de if y ariou ® levels of representational abstraction, but they 
should always have the three semiotic components so that the knowledge 
ke described and propagated in a manner similar to the 
actual cognitive process. This will ensure that the full range of 
engineering discourse, from the qualitative to the quantitative will 
be modeled by computer descriptions. M ' 11 

Furthermore, all three semiotic components are described by both a 
declarative and a procedural statement. The declarative statement 

covers »how” to "usJ “t?^ ^ dGSign ' and the P roc edural statement 


COMPONENTS OF AN ENGINEERING DEVICE 


LINGUISTIC 

ASPECTS 

ENGINEERING 

ASPECTS 

SEMANTIC 

FUNCTION 

SYNTACTIC 

STRUCTURE 

PRAGMATIC 

BEHAVIOR 

CONTEXT 


Figure 3 
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KNOWLEDGE STRUCTURE 


It is our contention that components of knowledge used in processes 
apparently as different as design and analysis are, in fact, the same. 
The description of each component and its processing vary as a 
function of the particular requirements of a problem situation, but 
the component itself stays the same. We propose that different design 
/ analysis reasoning procedures can, in fact, be represented as 
different refinements of the same higher-order semantic net. The 
various levels of detail required to solve problems correspond to 
various levels of representational abstraction. The same can be said 
for the representation of the facts in the domain of knowledge: 
Functional and structural hierarchies of the components of a complex 
system can be described at various levels of abstraction. 

We therefore propose the Multiple Layer Semantic Net (MLSN) [6] as the 
cognitive knowledge structure which unifies the representation of the 
various types of knowledge about facts and reasoning. The MLSN is 
conceptually a layered semantic net. The nets of each layer are 
isomorphous to one another in that they represent the same engineering 
concepts, but their descriptions of the concepts are made at different 
levels of abstraction. The descriptions are qualitative toward the top 
of the representation and quantitative toward the bottom. 

The rest of this presentation describes the cognitive techniques the 
MLSN should handle and points out the necessity to provide such a 
unified structure. 


QUALITATIVE 









MULTI-DISCIPLINARY ASPECTS OF THE DOMAIN OF KNOWLEDGE 


Most design problems require a combination of knowledge from different 
domains. For example, in the design of wood structures [7], wood 
science, wood engineering, and structural engineering are combined. In 
building design [8], it is architectural, structural, mechanical, and 
electrical engineering; in aerospace structures [9], aerodynamics, 
structural engineering, and mechanical engineering. In some design 
problems, the interaction among the various knowledge domains may be 
mostly sequential for the larger components of the process, whereas 
some sub-problems could be solved in parallel [10]. In all cases, a 
strong interaction exists among the different sources of knowledge, a 
fact which calls for new approaches such as simultaneous engineering 
and integrated activities. 

The complex structure being designed, e.g., a building, may be 
decomposed differently in each one of the knowledge domains and may 
have different function hierarchies in these domains. These various 
views of the same complex structure can be represented with 
corresponding hierarchies in the levels of the MLSN. The hierarchies 
of the different domains are interconnected by the appropriate 
semantic links, which account for the particular aspects of the 
context in which the complex structure is used within each discipline. 
An example is the relationship between the structural decomposition 
provided by an architect, which becomes the functional decomposition 
serving as the starting point for design by a structural engineer. 


COMPLEX STRUCTURE 



Figure 5 
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DECLARATIVE AMD PROCEDURAL DEFINITIONS 


Two different kinds of knowledge are used to perform a cognitive 
activity: Declarative knowledge and the procedural knowledge. 
Declarative knowledge consists of what we know about events, objects, 
and the relationships between them. Declarative knowledge is also 
referred to as propositional knowledge and can easily be represented 
by semantic networks [2, 11]. Procedural knowledge describes how to 
perform various activities and the dynamic process of how and why 
operations are performed upon the declarative knowledge. 

At a higher conceptual level, declarative and procedural descriptions 
are part of different knowledge processing skills. According to [12], 
we first form some declarative knowledge while learning a task; we 
then correct the declarative knowledge in the associative stage to 
form some procedural knowledge; in the autonomous stage, these 
procedures become highly automated. In familiar problems, experts use 
procedural knowledge in a relatively rapid and automatic fashion [13, 
14] and in a new and unusual situation they still have to rely on 
their declarative knowledge. 

Hence we propose that procedural knowledge is used for routine designs 

[15] , declarative knowledge is used for creative and innovative 
designs, and a combination of both is used for design by redesign 

[16] . 
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FROM QUALITATIVE TO QUANTITATIVE, COMPLETE OR PARTIAL ENVISIOMMENT 

The human mind can envision a complex system in its entirety or zero- 
in on one part of it. In doing so, it switches from higher levels of 
abstraction where the information tends to be more gualitative to 
lower levels of abstraction where it is more quantitative [6]. This is 
exemplified in the decompositional stage of design in which one 
critical component or sub-system is designed in more detail with the 
assumption that it will later fit with the rest. 

Just as the human mind shifts from qualitative to quantitative 
descriptions, so does the design process. A design at a deeper level 
of description defines one at a higher level by providing more detail 
about the components. This characterization corresponds to the 
cognitive process of definition (the reverse of abstraction) . It can 
also describe the reasons for having to define more precisely the 
concepts parametric in nature and includes the procedures to do so. 

Modeling the design knowledge in multiple layers is especially 
appropriate in routine design [15]: The structures being designed and 
their components stay fundamentally the same from one application to 
the next. Only the numerical values of the parameters change from one 
specialization to another. It is therefore not necessary to abstract 
toward the generalized conceptual structure, design after design. This 
process corresponds to moving from one level upward, then back down in 
the MLSN in a fundamentally qualitative-then-quantitative process. 


QUALITATIVE 
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/ STIFF-SYSTEM / 

/ HEAVY-LOADS / 

n / / 
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ft / / 

1 QUANTITATIVE 


LATERALLY STRESSED BRIDGE: 




27 ROUTINE DESIGNS WITH: 



SPAN < 50 FEET 
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SPAN = 24-30, 30-36, 36-42 FEET. 
ONE, TWO, OR THREE LANE-BRIDGE 
HS-15 , HS-20 , HS-25 LOADINGS 
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SIMULATION: DESIGN AMD EVALUATION 


The procedures of design and evaluation are dual of one another in the 
following sense: Design consists of creating the structure of a device 
that exhibits a specific and desired behavior or that is meant to 
serve an intended purpose. Evaluation, on the other hand, consists of 
analyzing the behavior of a device in an effort to understand what its 
structure must be for it to exhibit that behavior. Both design and 
evaluation processes use the same knowledge base of facts and 
relations; only the manipulations of the components vary between the 
processes, as will be shown later. 

Design and evaluation can be viewed as two refinements of the concept 
of simulation. Simulation is the attempt to make the composition of a 
system exhibit a certain behavior, and depends on the ability to 
create the system in the first place, whether it is a preliminary 
design alternative or a model of an existing system. 

Because of their duality and generalization to the same concept, it is 
logical to integrate a design and a risk assessment into the same 
program: The structure of a complex system is established to some 
degree of completeness during a preliminary design. That structure can 
then be investigated to evaluate the risk associated with a potential 
failure of some of the components of the structure. The decision to 
accept or reject the preliminary design alternative is then made based 
on the results of the risk analysis. 


DESIGN: BEHAVIOR OR FUNCTION =* STRUCTURE 
EVALUATION: STRUCTURE => BEHAVIOR AND/OR FUNCTION 

DESIGN AND EVALUATION ARE DUAL ON A SEMIOTIC BASIS 

CONSEQUENCES: . GENERALIZATION OF BOTH TO SIMULATION 

. INTEGRATION OF PRELIMINARY DESIGN AND 
RISK ASSESSMENT 



Figure 8 
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FUNDAMENTAL DESIGN STRATEGIES AND THEIR COGNITIVE EQUIVALENTS 

Every design involves four steps: Problem formulation, conceptual 
design, embodiment design, and final design. The first step of the 
conceptual design establishes the functional decomposition of a 
complex system and its components. This decomposition corresponds to 
the cognitive processes of 1 ) specialization of a concept into an 
instance, and 2) individuation of the concept into sub-components. 

° f the f on< r e P tual design is the design synthesis. This 
sembles some components into a more complex structural hierarchy 
which corresponds to the earlier functional decomposition. The 
corresponding cognitive process is the aggregation of concepts. 

Some basic design strategies applicable during the conceptual design 
are the routine design, design by redesign, innovative design, and 
creative design. Any combination of these can lead to even more 
complex strategies. 

Design by redesign first generalizes a concept to a higher-order 

an ? then s P eci alizes to another instance. Routine design 
first abstracts to a more qualitative model of the same structure and 

another more quantitative model. Both processes are 
sketched on the MLSN below. As already mentioned, procedural knowledqe 
is used in routine design, declarative knowledge in creative design 9 
and a combination of both in innovative design and design by redesign 
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Cognitive processes: 
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Design strategies and the MLSN: 

4 Routine design: path from si up to s, then down to s2 
. Design by redesign: path si to g' , then to s2 

Figure 9 
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DESIGN PROCEDURES AMD THE SEMIOTIC RELATIONSHIPS 

The reasoning procedures of the design problem solving process the 
knowledge among the components of a semiotic device, either by 
deriving one semiotic component from another inside one device, or by 
comparing similar components between two devices. There are six 
possible relationships among the three semiotic components of a 
device, all used either in design or analysis. 

The FUNCTION- to-STRUCTURE mapping (i.e., deriving the structure from 
the function) and the BEHAVIOR-to-STRUCTURE mapping take place in the 
design synthesis. They use teleological reasoning. The STRUCTURE-to- 
FUNCTION and the STRUCTURE-to-BEHAVIOR mappings are analysis 
processes. They use causal reasoning. 

Except for the mapping from structure to behavior, all mappings are of 
the type one to many. For example, several functions can be met by one 
structure, just as multiple structures could serve one function. A 

given structure can only generate one behavior at a time, with a fixed 
context . 


The FUNCTION- to-BEHAVIOR mapping can be part of the innovative design 
which consists of finding new applications to an existing device. This 
mapping can be one to many. Finally, the BEHAVIOR-to-FUNCTION mapping 
corresponds to a qualitative analysis process and is a one to one 
mapping if considered in one context. 



design 


Figure 10 
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DESIGN STRATEGIES AND THE SEMIOTIC RELATIONSHIPS 

The design process is based at the fundamental level on causal and 
teleological reasoning. Causal reasoning processes "what it is" in 
order to derive "what it does". It is applied, for example, in a 
backward chaining manner in the FUNCTION-TO-BEHAVIOR mapping of an 
innovative design where a new usage is identified for a device. 
Teleological reasoning, by contrast, processes "what it is for" to 
derive "what it should be". It is applied, for example, in the 
traditional derivation of the STRUCTURE from the FUNCTION. 

At a higher level, some design strategies are Design by analogy, 
which compares corresponding components of different devices; design 
by constraint satisfaction, which builds up information requirements 
from the context for the function and structure of a device; and 
design by analysis, as in the innovative design process mentioned 
above. In case-based designs as in design by analogy, all 
transformations could be used [ 17 ] . 

Even higher order design strategies still manipulate the semiotic 
components. Routine design involves transformations of a structure 
from one instance into another one. Design by redesign involves 
iterations on the transformations between the function and/or the 
behavior and the structure. In all multidisciplinary designs, 
structures of one domain are functions for another. Through the design 
process, the structures of the second domain finish completing the 
description of the initial structures. 


REASONING, PROCESSES, AND TYPES OF DESIGN 
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SUMMARY 


T tu 5 ln t U ^° matlng the desi 9 n and evaluation of complex systems 
led to the formulation of a cognitive knowledge structure developed 
to facilitate the acquisition and representation of knowledge at 
multiple levels of abstraction. 

The knowledge structure, the multiple layer semantic nets (MLSN) 
consists of isomorphous semantic nets describing the relationships 
among concepts viewed as semiotic paradigms. The components of the 
semiotic paradigms (structure, function, behavior and context) are 
described from qualitative levels to quantitative levels by both 
declarative and procedural descriptions. 

T?* des f ribed . h ere in the perpective of the design process 

and the design strategies it should handle. It is also applied in 
another component of the research to investigate and develop 

systems 1163 ' baSed ° n qualitative reasoning, to evaluate complex 


The MLSN is now used to guide the development of 
which will perform both the design and the risk 
structural systems. 


a computer 
assessment 


program 
for complex 
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